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ABSTRACT 


This  thesis  investigates  the  application  of  Wiener  filtering  and  wavelet 
techniques  for  the  removal  of  noise  from  underwater  acoustic  signals.  Both  FIR  and 
HR  Wiener  filters  are  applied  in  separate  methods  which  involve  the  filtering  of 
Wavelet  coefficients  which  have  been  produced  through  a  discrete  wavelet 
decomposition  of  the  acoustic  signal.  The  effectiveness  of  the  noise  removal 
methods  is  evaluated  by  applying  them  to  simulated  data.  The  combined  Wiener 
wavelet  filtering  methods  are  compared  to  more  traditional  denoising  techniques 
which  include  Wiener  filtering  and  wavelet  thresholding  methods. 
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I.  INTRODUCTION 


Accurate  analysis  of  ocean  acoustic  signals  has  proven  to  be  a  difficult  task  due  to 
unwanted  noise  which  is  usually  present.  Significant  effort  has  been  directed  to  the  problem 
of  noise  removal  from  underwater  acoustic  signals.  In  the  1940s,  Norbert  Wiener  conducted 
extensive  research  in  the  area  of  linear  minimum  mean-square  error  (MMSE)  filtering. 
Provided  the  spectral  characteristics  of  an  additive  combination  of  signal  and  noise  are 
known,  the  Wiener  filter  will  operate  in  a  linear  manner  to  yield  the  best  separation  of  the 
signal  from  the  noise.  Here,  “best”  means  minimum  mean-square  error.  The  Wiener  filter 
is  today  a  well  tested  method  of  noise  reduction. 

Wavelet  techniques  are,  conversely,  a  more  modem  approach  to  noise  reduction, 
although  their  origin  lies  in  the  timeless  methods  developed  by  Fourier.  Wavelet  analysis 
decomposes  the  signal  into  a  family  of  basis  functions  and  provides  two  significant 
advantages  over  the  traditional  Fourier  analysis.  First,  in  wavelet  analysis  there  is  a  wide 
variety  of  basis  functions  to  choose  from,  and  second,  a  multi-resolution  capability  is 
provided  in  the  time-frequency  domain  which  is  critical  to  the  identification  and  elimination 
of  noise  in  a  non-stationary  environment.  Wavelet  techniques  have  proven  to  be  a  viable  tool 
for  the  denoising  of  acoustic  signals. 

Recently,  a  question  has  been  posed  concerning  the  combination  of  these  two 
invaluable  techniques.  The  approach  presented  applies  a  discrete  wavelet  transform  to  the 
noisy  signal.  Next,  rather  than  thresholding  the  wavelet  coefficients,  a  Wiener  filter  is 
applied  to  separate  the  noise  from  the  signal.  While  this  approach  may  intuitively  seem 
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reasonable,  there  are  problems  associated  with  aliasing  and  perfect  reconstruction  of  the 
denoised  signal  which  must  be  considered.  Given  that  the  aliasing  problem  can  be  dealt 
with,  the  results  still  may  not  be  superior  to  either  Wiener  filtering  alone  or  wavelet-based 
techniques  alone  but  the  possibilities  are  certainly  worthy  of  investigation. 

In  this  thesis,  an  approach  is  presented  which  combines  the  Wiener  filter,  in  both  the 
causal  finite  impulse  response  (FIR)  and  the  non-causal  infinite  impulse  response  (HR) 
forms,  with  wavelet-based  techniques.  The  various  noise  removal  methods  are  applied  to 
simulated  data  which  offers  the  advantage  of  providing  ground  truth  to  the  analysis. 
Additionally,  the  signal-to-noise  ratio  (SNR)  level  can  be  easily  modified  to  evaluate  the 
effectiveness  of  the  denoising  algorithms.  The  combined  Wiener  wavelet  based  denoising 
methods  demonstrate  promise  in  the  recovery  of  ocean  acoustic  signals  from  the  noisy  ocean 
environment  in  which  ships  operate. 

The  problem  of  denoising  underwater  acoustic  signals  is  addressed  in  nine  additional 
chapters  summarized  as  follows.  Chapter  II  presents  background  information  on  the 
characteristics  of  ocean  acoustic  noise .  In  Chapter  HI  the  theory  and  application  of  Wiener 
filtering  methods  is  discussed.  Chapter  IV  presents  the  theoretical  development  of  wavelet 
analysis  and  and  mutirate  systems.  Standard  wavelet  threshold  based  denoising  techniques 
are  presented  in  Chapter  V.  FDR.  Wiener  wavelet  noise  removal  methods  are  developed  in 
Chapter  VI  followed  by  DOR  Wiener  wavelet  methods  in  Chapter  VII.  A  comparison  of  the 
various  denoising  methods  is  presented  in  Chapter  VIH.  Wavelet  packet  denoising 
techniques  are  applied  to  the  problem  of  high  frequency  signals  of  short  duration  in  Chapter 
IX.  Finally,  conclusions  are  presented  in  Chapter  X. 


2 


II.  OCEAN  ACOUSTIC  NOISE 


Though  significant  resources  have  been  devoted  to  non-acoustic  detection  methods 

in  the  undersea  warfare  environment,  sonar  remains  the  primary  and  most  reliable  method 

of  detection.  However,  as  modem  technology  continues  to  enable  significant  reduction  in 

radiated  source  levels  the  problem  of  early  detection  and  classification  of  underwater 

acoustic  signals  gains  greater  significance.  A  possible  solution  lies  in  the  area  of  improved 

digital  signal  processing  and  filtering  methods  which  would  allow  detection  and 

classification  of  signals  at  lower  signal-to-noise  ratios.  The  passive  sonar  equation  states  that 

the  source  level  of  the  target  minus  the  loss  due  to  propagation  through  the  medium,  minus 

the  sum  of  all  interfering  noises  plus  improvement  by  the  spatial  processing  gain  of  the 

receiver,  must  be  greater  than  the  detection  threshold  for  a  sonar  system  to  sense  the  presence 

of  a  target  with  a  fifty  percent  probability  of  detection  [1] 

SL  -  TL  >  NL  -  DI  +  DT,  (2.1) 

where:  SL  =  Source  Level  of  the  target  being  detected  passively 

TL  =  Transmission  Loss  as  the  signal  propagates  to  the  detector 

NL  =  Noise  Spectrum  Level  of  the  ambient  noise  and  self  noise  in  the  ocean 

DI  =  Directivity  of  the  detecting  system 

DT  =  Detection  threshold  for  50%  probability  of  detection 

and  all  terms  are  in  dB  referenced  to  lp  Pa. 

A  reduction  in  the  noise  spectrum  level  would  certainly  result  in  an  improved 

detection  probability.  The  noise  spectrum  level  includes  both  self  noise  and  ambient  noise. 
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A.  NOISE  SOURCES 

1.  Ambient  Noise 

Urick  identifies  ambient  noise  as  that  part  of  the  total  noise  background  observed 
with  a  nondirectional  hydrophone  which  is  not  due  to  that  hydrophone  and  its  manner  of 
mounting  called  “self-noise”,  or  to  some  identifiable  localized  noise  source  [2].  Ambient 
noise  is  produced  by  a  variety  of  sources  and  may  be  found  in  the  frequency  range  from  1 
Hz  to  100  KHz.  The  noise  frequency  range  is  typically  divided  into  five  frequency  bands  of 
interest,  which  are  discussed  next. 

a.  Ultra-Low  Band  ( <1  Hz) 

Though  little  is  known  about  the  exact  contributors  at  the  lower  end  of  the 
spectrum,  it  is  surmised  that  these  sources  include  tides  and  hydrostatic  effects  of  waves, 
seismic  disturbances,  and  oceanic  turbulence  [3].  Tides  and  waves  cause  hydrostatic 
pressure  changes  resulting  in  a  low  frequency,  high  amplitude  contribution  to  the  ambient 
noise  spectrum.  The  constant  seismic  activity  measured  on  land  extends  into  the  ocean 
environment  causing  low  frequency,  high  amplitude  contributions  which  add  to  those 
produced  by  tides  and  waves.  Oceanic  turbulence  gives  rise  to  varying  dynamic  pressures 
which  are  detected  by  pressure  sensitive  hydrophones. 

b.  Infrasonic  Band  (1  to  20  Hz) 

This  band  has  gained  importance  with  the  emergence  of  improved  low 
frequency  narrowband  passive  sonar  systems.  It  contains  the  strong  blade-rate  fundamental 
frequency  of  propeller-driven  vessels  and  its  accompanying  harmonics.  A  steep  negative 
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spectral  slope  of  10  dB  per  octave  is  common  in  the  region  from  1  to  5  Hz.  This  slope  goes 
positive  from  5  to  20  Hz  as  shipping  noise  begins  to  become  a  more  significant  factor.  In  the 
absence  of  ship  traffic  this  region  continues  to  fall  off  and  is  more  affected  by  wind  speed. 

c.  Low  Sonic  Band  (20  to  200  Hz) 

Studies  have  shown  that  distant  ship  traffic  is  the  dominant  source  of  noise 
at  100  Hz  and  has  a  significant  effect  in  the  low  sonic  band  [3].  In  areas  of  low  shipping 
intensity,  wind  speed  continues  to  be  the  major  factor  just  as  it  is  in  the  infrasonic  and  high 
sonic  bands.  Thus,  an  area  of  heavy  shipping  such  as  the  North  Atlantic,  where  on  average 
1,100  ships  are  underway,  will  see  a  much  greater  effect  than  less  traveled  areas  such  as  the 
South  Pacific. 

d.  High  Sonic  Band  (200  to  50,000  Hz) 

The  well  known  acoustician  V.O.  Knudsen  conducted  extensive  studies  in  this 
band  during  World  War  n.  In  these  studies,  he  was  able  to  correlate  noise  with  wind  speed 
in  the  frequency  range  500  Hz  to  25  KHz.  His  results,  published  in  1948,  are  best 
summarized  by  the  curves  shown  in  Figure  2. 1 .  These  curves  show  a  clear  relationship 
between  wind  speed  (or  sea  state  which  isn’t  measured  as  accurately)  and  spectrum  levels. 
Subsequent  studies  have  shown  the  spectrum  to  be  flatter  in  the  range  200  to  800  Hz  but 
have  generally  confirmed  Knudsen ’s  results  [4].  Crouch  and  Burt  [5]  have  developed  an 
expression  to  model  the  noise  spectrum  level  in  dB  which  is  given  as 

NL(f)  =  B(f)  +  20  n  log10  V,  (2.2) 

where  NL  is  the  noise  spectrum  level  in  dB  referenced  to  1  pPa  at  frequency/,  B(f)  is  the 
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noise  level  at  a  wind  speed  of  1  knot  at  a  particular  frequency,  n  is  an  empirical  coefficient, 


Figure  2.1:  Average  Ambient  Noise  Spectra.  From  Ref.  [6]. 

and  V  is  the  wind  speed  in  knots.  For  n  =  1  the  noise  level  increases  as  20  log10  V,  and  the 
noise  intensity  will  increase  as  the  square  of  the  wind  speed. 
e.  Ultrasonic  Band  (>  50,000  Hz) 

At  frequencies  above  50,000  Hz  thermal  noise  becomes  the  predominant 
contributor  to  the  noise  background.  In  1952  Mellen  [7]  showed  theoretically  that  the 
thermal  noise  of  the  molecules  of  the  sea  affects  hydrophones  and  places  a  limit  on  their 
sensitivity  at  high  frequencies.  Based  on  principles  of  classical  statistical  mechanics  he 
formulated  the  following  expression  for  the  spectrum  level  NL  in  dB  referenced  to  1  pPa  of 
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the  thermal  noise  at  frequency/  in  kHz  given  by 

NL(f)  =  -15  +  20  logj0/  (2.3) 

Though  measurements  have  been  recorded  in  this  band  they  have  not 
conclusively  substantiated  the  above  expression  due  to  excessive  shipping  noise  in  nearby 
ports.  No  measurements  in  this  band  in  deep,  quiet  open  ocean  water  appear  to  have  been 
made  until  now. 

2.  Self  Noise 

Self  noise  includes  all  noise  created  by  the  receiving  platform  and  usually  falls  into 
one  of  two  categories  which  include  equipment  noise  and  platform  noise.  Equipment  noise 
includes  electronic  or  thermal  noise  produced  within  the  sonar  electronic  system.  Platform 
noise  is  produced  from  the  same  sources  as  radiated  noise  except  that  the  source  of  platform 
noise  is  the  receiving  platform  vice  the  source  or  target  platform.  These  sources  include 
machinery  noise,  hydrodynamic  noise,  propeller  noise  and  transients.  Platform  noise  reaches 
the  receiving  transducer  by  a  variety  of  methods  including  vibration  via  an  all-hull  path,  all¬ 
water  direct  path,  all-water  back  scattered  path  from  volume  scatterers,  and  all-water  bottom 
reflected  path  [2].  Machinery  noise  occurs  principally  as  low  frequency  tonals  which  are 
relatively  independent  of  speed.  Hydrodynamic  noise  which  includes  all  sources  of  noise 
resulting  from  the  flow  of  water  past  the  hydrophone  and  its  support  and  the  outer  hull 
structure  of  the  vessel,  becomes  more  significant  as  speed  increases.  At  high  speeds 
propeller  noise  becomes  the  dominant  contributor  to  self  noise. 
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B.  TRANSMISSION  LOSS  AND  WATER  MASS  CHARACTERISTICS 

Properties  of  the  water  mass  such  as  temperature,  salinity,  and  density  affect  the 
transmission  path  of  sound  in  water  and  therefore  alter  the  signal  received  at  the  hydrophone. 
Additionally,  the  depth  of  water  and  bottom  structure  influence  the  path  traveled  and  could 
result  in  multiple  transmission  paths  between  the  source  and  receiver.  The  affects  of 
absorption  and  attenuation  cause  the  ocean  to  filter  out  the  high  frequency  spectrum  while 
passing  the  low  frequency  spectrum.  Thus,  ocean  acoustic  noise  appears  to  occur  more 
predominantly  in  the  lower  frequency  region. 

C.  NOISE  MODEL 

Modeling  ocean  acoustic  noise  as  a  stationary  white  Gaussian  random  process 
substantially  simplifies  the  analysis  and  study  of  denoising.  However,  the  variety  of  sources 
which  contribute  to  self  and  ambient  noise  and  the  acoustic  transmission  characteristics  of 
the  ocean  result  in  a  noise  contribution  which  is  colored  vice  white  in  nature.  A  common 
method  around  this  obstacle  is  to  pre-whiten  the  acoustic  signal. 

Stationarity  is  another  assumption  which  must  be  applied  carefully.  Certainly  the 
ocean  environment  is  one  in  which  sources  of  acoustic  signals  change  regularly.  Moreover, 
the  water  mass  properties  are  in  a  constant  state  of  flux.  This  hardly  meets  the  criterion  for 
stationarity  at  first  glance.  However,  on  a  time  scale  of  several  seconds,  the  ocean 
environment  can  indeed  be  considered  stationary. 

Finally,  evaluation  of  ocean  acoustic  noise  reveals  that  the  assumption  of  a  Gaussian 
random  process  appears  to  hold  true  in  most  cases.  An  artificial  computer  generated  white 
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Gaussian  noise  sample  and  its  power  spectral  density  are  pictured  in  Figure  2.2.  This 
assumption  can  be  proved  by  comparing  the  histogram  of  a  noise  sample  to  the  Gaussian 
probability  density  function  with  appropriate  sample  mean  and  variance  or  by  checking  a 
normal  probability  distribution  plot  of  the  noise  sample  as  seen  in  Figure  2.3.  Barsanti’s 
results  have  verified  this  assumption  with  many  actual  ocean  acoustic  signals  [8],  Frack 
discusses  several  more  involved  methods  of  verifying  this  critical  modeling  assumption  [9]. 

All  noisy  signals  in  this  thesis  are  generated  by  adding  white  Gaussian  noise  with 
zero  mean  (as  shown  in  Figure  2.2)  to  the  noise-free  signal.  When  analyzing  actual  ocean 
acoustic  signals  a  noise  sample  must  be  available  to  determine  the  statistical  properties 
necessary  to  produce  optimal  filtering.  A  noise  sample  with  statistical  properties  which 
accurately  reflect  the  noise  embedded  in  the  noisy  signal  produces  a  more  effective  filter  and 
a  more  accurate  denoised  signal.  Frack  discusses  methods  which  can  be  used  to  model  noise 
for  filtering  purposes  when  adequate  noise  samples  are  not  available  [9]. 
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Normal  Probability  Plot  of  Noise 


Figure  2.3:  Histogram  and  Normal  Probability  Plot  of  Artificial  White  Gaussian  Noise 
Sample.  The  superimposed  normal  PDF  has  been  scaled  to  match  the  histogram  for  1024 
samples. 
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III.  WIENER  FILTERING 


Wiener  filtering  is  one  of  the  earliest  methods  used  to  separate  noise  from  a  desired 
signal.  Developed  by  Norbert  Wiener  in  the  1940's,  this  form  of  optimal  filtering  is  used  to 
denoise  ocean  acoustic  signals  with  the  objective  of  improving  signal  detection  and 
classification.  The  FIR  Wiener  filter  and  its  application  to  denoising  is  developed  in  this 
chapter. 

A.  MODEL  DESCRIPTION 

The  optimum  linear  discrete-time  estimation  model  is  illustrated  below  in  Figure  3. 1 . 
The  filter  input  signal  x(ri),  consists  of  a  signal  s(n)  and  an  additive  noise  w(n).  The 
estimator  is  constrained  to  be  a  linear  filter  with  impulse  response  h(n),  which  is  designed 
to  remove  the  noise  while  preserving  the  characteristics  of  the  signal.  The  filter  impulse 
response  is  obtained  by  minimizing  the  estimation  error  e(n ),  defined  as  the  difference 
between  the  filter  output  and  the  desired  signal  d(ri),  which  is  taken  as  the  original  signal 
s(n) ,  for  filtering  applications. 
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A  time- varying  filter  can  be  derived  for  nonstationary  signals  but  the  complications 
are  significant  [10].  A  more  tractable  approach  involves  segmenting  the  input  signal  into 
blocks  which  can  be  considered  stationary  for  the  short  duration  of  a  few  milliseconds  to  a 
few  seconds.  In  developing  the  filter  for  this  model,  the  input  signal  and  noise  spectral 
characteristics  or  equivalently,  auto-  and  cross-correlation  functions  must  be  known. 
Although  typically  not  known  a  priori,  this  information  can  be  estimated  from  the  data 
assuming  that  the  noise  can  be  isolated  sometime  during  the  observation  interval. 

B.  FIR  WIENER  FILTER  DEVELOPMENT 

Given  the  spectral  characteristics  of  an  additive  combination  of  signal  and  noise, 
Wiener  proposed  a  scheme  which  best  separates  signal  from  the  noise  [11].  His  procedure 
involves  solving  for  the  “best”  filter  coefficients  by  minimizing  the  difference  between  the 
desired  ideal  noise-free  signal  and  the  filter  output  in  a  mean  square  error  sense.  The  mean 
square  error  is  most  often  used  due  to  its  mathematical  simplicity,  as  it  leads  to  a  second 
order  cost  function  with  a  unique  minimum  obtained  by  using  optimal  filter  coefficients  [10]. 
The  FIR  Wiener  filter  is  constrained  to  length  P  with  coefficients  hk  (0  <  k  <  P-1). 


Its  output  depends  on  the  finite  input  data  record x(n)  =  [x(n),  x(n-l), . x(n-P+l)]T  and  may 

be  expressed  as  the  convolutional  sum: 

p-i 

s(n )  =  y(n)  =  ^h  *(k)x(n  -k).  (3.1) 

k= 0 


The  mean  square  value  of  the  error,  e(ri),  between  the  desired  output  s(ri),  and  the  filter 
output  s(n) ,  can  then  be  expressed  as: 

p- 1 

ol=E{\e(n)\2}  =£{ls(n)-£  h  \k)x(n-k) I2}.  (3.2) 

k=0 
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The  principle  of  orthogonality  states  that  the  optimal  filter  coefficients  h(ri),  for  n  =  0,1,..., 
P-1,  minimize  the  mean  square  error  if  chosen  such  that  E{x(n-i )  e\n)}  =  0,  for  i  =  0,1,..., 
P-1,  that  is  if  the  error  is  orthogonal  to  the  observations  [12].  The  minimum  mean  square 
error  can  then  be  given  by  o2e  =  E{s(n)  e\n)}.  Now  applying  the  orthogonality  principle 
results  in  the  following  set  of  equations: 


E{x{n-i)e  *(n)}  =E\ 


'  p- i 

*  x(n-i ) 

s  h(k)x  *(n-k )  " 

*=o  J 

r  =  0,  for  i  =  0,l,. ..,P-1.  (3.3) 


Equation  (3.3)  can  also  be  expressed  in  terms  of  the  auto-correlation  function  of  the 

observations  Rx{k)  and  the  cross-correlation  function  between  the  signal  and  the  observation 

sequence  Rsx(i)  which  leads  to: 
p- i 

'£Rx(k-i)h(k)  =  Rsx*(i),  for  i  =  0,1,...,P-1.  (3.4) 

k=0 

Assuming  the  signal  and  noise  are  independent  and  have  a  zero  mean,  it  can  be  shown  that: 

RJk)=Rs(k)  (3.5) 


and  Rx(k)=Rs(k)+Rw(k), 

where  Rs(k)  and  Rw(k)  are  the  signal  and  noise  correlation  functions  respectively. 

The  set  of  discrete  Wiener-Hopf  equations  for  the  causal  stationary  filter  can  be 
expressed  in  matrix  form  as 


RXh  =  fsx’ 

(3.6) 

where 

Rx  =  E{x(n)x'T(n)} 

(3.7) 

and 

rsx  =  E{s(n)x(n)} 

(3.8) 

and  represents  the  reversal  operator  applied  to  a  vector  [12].  Equation  (3.6)  may  be 
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solved  for  the  filter  coefficients  using  matrix  methods.  The  minimum  mean  square  error  is 
found  using  the  second  part  of  the  orthogonality  theorem: 


Og  =  E{s(ri)e*{ri)}  =  Zsj  s(n)|^  s  h  \k)x{n-k) 

or  again  in  terms  of  the  correlation  function: 

o\  =  R(0)-Y,h*(k)Rsx(k). 


(3.9) 


(3.10) 


k= 0 

The  Wiener  filter  is  now  applied  to  a  synthetically  generated  noisy  sinusoidal  signal 
of  varying  frequency  with  a  SNR  of  0  dB.  The  original  signal  and  noisy  signal,  shown  in 
Figure  3.2(a),  are  compared  to  the  original  signal  and  denoised  signal  in  Figure  3.2(b).  The 
standard  measure  used  to  compare  denoised  signals  to  the  noise  free  original  signal  is  a 


modified  Mean  Squared  Error  ( MSE ),  defined  as: 


MSE 


£ 


n- 1 


s{n)  y(n) 
norm(s )  normiy) 


2 


(3.11) 


where  s(n)  is  the  noise  free  original  signal,  y(n)  is  the  denoised  output,  and  N  is  the  length 
of  the  signal.  The  signals  are  energy  normalized  by  dividing  by  their  norms.  This  represents 
a  denoised  signal  amplitude  gain  which  is  applied  to  account  for  losses  incurred  during  the 
filtering  process.  This  normalized  MSE  will  be  used  throughout  the  thesis  as  a  measure  of 
performance.  It  is  not  normalized  by  the  signal  length.  However,  a  standard  signal  length 
of  1024  or  16384  points  is  used  exclusively  for  all  denoising  trials.  The  Wiener  filtering 
results  for  various  filter  orders  is  depicted  in  Figures  3.3  and  3.4.  A  single-frequency  noisy 
sinusoidal  signal  ( 0  dB  SNR)  is  denoised  through  ten  trials  and  an  average  normalized  MSE 
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is  computed.  This  is  repeated  for  test  signals  with  normalized  frequencies  which  range  from 
0.01  to  0.49  in  steps  of  0.01  and  the  results  are  displayed.  The  noise  sample  supplied  to  the 
Wiener  filter  is  an  independent  noise  sample  in  Figure  3.3  and  the  original  noise  sample  used 
to  create  the  noisy  sinusoidal  signal  in  Figure  3.4.  Although  the  original  noise  sample  is  not 
practically  available,  it  is  shown  for  comparison  purposes.  A  table  comparing  the  MSE 
values  averaged  across  the  spectrum  is  shown  following  the  Figures.  The  MSE  performance 
improves  for  higher  filter  orders.  Additionally,  the  MSE  value  is  rather  flat  across  the 
spectrum  for  a  filter  orders  higher  than  four.  These  results  are  as  expected  for  the  Wiener 
filter  and  provide  a  benchmark  standard  of  performance  against  which  other  denoising 
methods  can  be  compared. 
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Amplitude  Amplitude 


Original  Signal  and  Noisy  Signal:  SNR  =  OdB,  Normalized  Frag  =  0.1 


Xime  (sample  number) 


Original  Signal  and  Wiener  Filtered  Signal:  Filter  Order  =12 


Xime  (sample  number) 

Figure  3.2:  Wiener  Filter  (a)  Original  (dashed)  and  Noisy  Sinusoidal  Signal  at  normalized 
frequency=0.1,  sampling  frequency=l.  (b)  Original  (dashed)  and  filtered  signal. 
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Mean  Square  Error 


Wienr  Filter  Results:  SNR  =  OdB 


Independent  noise  sample.  Ten  trial  average. 
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Wiener  F"  lit  or  Results:  SNR  —  OdB 


Figure  3.4:  Mean-square  error  comparison  for  Wiener  filters  of  various  order.  Actual 
noise  sample.  Ten  trial  average. 


Wiener 

Filter  Order 

Actual  Noise 
Average  MSE 

Independent  Noise 
Average  MSE 

4 

0.3400 

0.3423 

8 

0.2078 

0.2167 

12 

0.1500 

0.1647 

20 

0.0977 

0.1275 

30 

0.0711 

0.1200 

18 


IV.  WAVELET  ANALYSIS  AND  MULTIRATE  SYSTEMS 


Wavelet  analysis  has  found  a  myriad  of  applications  in  areas  ranging  from  financial 
research  to  engineering  analysis.  It  is  particularly  well  suited  for  signal  processing 
applications  such  as  image  compression,  sound  enhancement,  and  statistical  analysis  [13]. 
This  chapter  presents  the  underlying  theory  of  wavelet  analysis  including  both  the 
mathematical  basis,  which  stems  from  Fourier  analysis,  and  the  multirate  implementation  in 
filter  banks. 

A.  FOURIER  ANALYSIS 

Fourier  and  wavelet  analysis  are  similar  in  that  the  goal  is  to  express  an  observed  real 
time-series  signal,  x(t),  as  a  linear  decomposition  of  the  form 

x(t)  =  (4-1) 

k 

where  k  is  an  integer  index  for  the  finite  or  infinite  sum,  ak  are  real-valued  expansion 
coefficients,  and  ek(t)  are  a  set  of  functions  called  the  expansion  set  [14].  If  an  appropriate 
expansion  set  is  chosen,  the  resulting  expansion  coefficients  will  provide  useful  information 
about  the  signal  allowing  better  analysis  and  processing.  The  expansion  set  is  called  a  basis 
for  the  signal  space  if  every  signal  in  that  space  can  be  expressed  as  a  unique  set  of  linearly 
independent  functions.  Further,  if  the  basis  functions  are  orthogonal  so  that  their  inner 
product  is  zero: 

<€//),  €k(t))  =  f  eft)ek{t)dt  =  0,  j*k  (4.2) 

then  the  coefficients  can  be  calculated  by  the  inner  product 
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ak  =  (x(t),ek(t))  =  Jx(t)e*k(t)dt. 


(4.3) 


In  addition  to  simplifying  the  computation  of  the  expansion  coefficients,  an  expansion 
utilizing  an  orthogonal  basis  function  results  in  a  sparse  representation  where  most  of  the 
coefficients  are  zero  or  near  zero.  This  has  great  utility  in  the  nonlinear  noise  reduction 
methods  developed  in  later  chapters. 

A  well  known  transform  which  is  used  to  extract  spectral  information  from  time- 
series  signals  is  the  Fourier  transform, 

oo 

X(f)  =  I  x(t)e  -2Jn/tdt.  (4.4) 

-oo 

which  uses  orthogonal  sinusoidal  basis  functions  contained  in  the  complex  exponential 
function.  Since  the  Fourier  transform  is  integrated  over  all  time,  it  is  only  suitable  for 
stationary  signals  which  have  no  time  varying  frequencies.  However,  most  real-world  signals 
are  non-stationary,  and  therefore,  a  time-varying  spectral  representation  is  essential  in 
conceptualizing  the  temporal  relationship  between  the  spectral  components. 

The  short-time  Fourier  transform  (STFT)  is  the  key  to  adapting  the  Fourier  transform 
to  time-frequency  analysis.  Let  x(t)  be  the  nonstationary  signal  of  interest.  A  temporal 
window,  w(t- x),  is  applied  to  the  signal.  The  window,  which  is  centered  at  x,  segments  the 
signal  into  short  intervals  which  can  be  considered  stationary.  The  traditional  Fourier 
transform  is  then  applied  to  the  windowed  signal.  Thus,  the  STFT  may  be  expressed  as 
follows: 

oo 

STFT(x,f )  =  f  xit)  w(t-x)  exp(  -j2nfi)dt.  (4.5) 

— oo 

As  a  result,  the  STFT  maps  a  one-dimensional  time  series  into  a  two-dimensional 
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time-frequency  function  STFT(x,f).  Equation  (4.5)  can  be  considered  the  inner  (scalar) 
product  of  the  signal  jc(t)  with  a  two  parameter  family  of  basis  functions  w(t- 1)  exp(-j27t ft) 
for  varying  x  and  /  [15].  The  squared  modulus  of  the  short-term  Fourier  transform  of  a 
signal  is  known  as  a  spectrogram  and  is  expressed  as  SPEC(x/)  =  I  STFT(x/)  1 2.  The 
spectrogram  is  a  measure  of  the  signal  energy  in  the  time-frequency  plane  and  provides 
powerful  insight  when  analyzing  nonstationary  signals.  A  common  band-pass  filter 
implementation  of  the  STFT  is  shown  in  Figure  4.1.  In  this  model,  the  frequency  /  might 
be  considered  the  midband  frequency  of  the  bandpass  filter.  This  model  results  in  a  filter 
bank  setup  with  a  set  of  analysis  filters  that  perform  a  STFT  for  each  specific  value  of 
frequency  /  [15]. 


Band-pass  Filter 


Figure  4.1:  Band-pass  Filter  Implementation  of  STFT. 

In  terms  of  time  and  frequency  resolution  the  STFT  must  satisfy  the  time-bandwidth 
limitations  of  the  uncertainty  principle,  namely: 


where  At  is  the  time  resolution  and  Af  is  the  frequency  resolution  or  bandwidth.  This 
equality  can  be  met  by  using  a  Gaussian  window  [16].  However,  once  the  window  is 
specified  a  constant  bandwidth  filter  bank  results  with  fixed  time-frequency  resolution  over 
the  entire  time-frequency  plane.  The  time-frequency  plane  has  a  uniform  tiling  as 
shown  in  Figure  4.2.  Although  the  Gaussian  window  minimizes  the  time-bandwidth  product, 
the  STFT  still  doesn’t  provide  sufficient  time  resolution  without  a  significant  degradation 
in  frequency  resolution.  A  narrow  or  compactly  supported  window  function  provides  better 
time  resolution  at  the  expense  of  poorer  frequency  resolution.  Similarly,  a  wider  window 
function  results  in  better  frequency  resolution  but  poorer  time  resolution.  An  alternative 
approach  to  the  STFT  is  required,  one  which  analyzes  the  signal  at  different  frequencies  with 
different  resolutions. 

B.  THE  CONTINUOUS  WAVELET  TRANSFORM 

The  wavelet  transform  is  a  multi-resolution  analysis  (MRA)  method  designed  to 


STFT  Time-Frequency  Display 


Figure  4.2:  STFT  Time-Frequency  Display. 
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provide  good  time  resolution  and  poor  frequency  resolution  at  high  frequencies  and  good 
frequency  resolution  and  poor  time  resolution  at  low  frequencies.  Fortunately,  many 
practical  ocean  acoustic  signals  are  composed  of  high  frequency  components  of  short 
duration  plus  low  frequency  components  of  long  duration,  and  thus  are  well  suited  to  this 
form  of  analysis  [17].  In  terms  of  the  filter  bank  model  described  earlier,  the  time  resolution 
must  increase  as  the  central  frequency  of  the  analysis  filters  increases.  The  filterbank  is  then 
composed  of  band-pass  filters  with  constant  relative  bandwidth  otherwise  known  as  constant- 
Q  filters  where  the  constant  relative  bandwidth  is  defined  as: 

=  constant  =  Q.  (4.7) 

A  clearer  understanding  may  be  obtained  by  comparing  the  division  of  the  frequency 
domain,  as  shown  in  Figure  4.3  below.  The  STFT  analysis  filters  in  Figure  4.3(a)  are 
regularly  spaced  over  the  frequency  axis  while  WT  multi-resolution  analysis  filters  in  Figure 


STFT  with  Constant  Bandwidth  Filters 


Figure  4.3(a):  STFT  decomposition  with  Constant  Bandwidth  Filters. 
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WTF  with  Constant-Q  Factor  Filters 


Frequency 


Figure  4.3(b):  WT  decomposition  with  Constant-Q  Factor  Filters. 

4.3(b)  are  regularly  spread  in  a  logarithmic  scale.  The  improved  frequency  resolution  at  low 

frequencies  and  better  time  resolution  at  higher  frequencies  is  clearly  illustrated  for  the  WT 

case. 

The  continuous  wavelet  transform  (CWT)  is  a  powerful  multi-resolution  analysis 
method  which  has  become  a  widely  used  alternative  to  the  STFT.  It  may  be  viewed  as  a 
signal  decomposition  onto  a  set  of  basis  functions  where  the  basis  functions  are  called 
wavelets  [15].  The  CWT  has  some  similarities  to  the  STFT  in  that  a  segmented  time-domain 
signal  is  transformed  by  multiplying  it  by  a  wavelet  function  which  is  similar  to  a  windowed 
function  in  the  STFT.  However,  the  basis  functions  for  the  two  transforms  are  quite 


Figure  4.4:  Basis  Functions  -  Sine  Wave  and  Daubechies-8 
Wavelet. 
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different.  Whereas,  the  Fourier  transform  breaks  the  signal  into  a  series  of  sine  waves  of 
different  frequencies,  the  wavelet  transform  breaks  the  signal  into  its  “wavelets”,  scaled  and 
shifted  versions  of  the  “mother  wavelet”. 

The  two  basis  functions  are  illustrated  in  Figure  4.4,  where  a  sine  wave  and  a 
Daubechies-8  wavelet  are  compared.  The  smooth  sine  wave  of  finite  length  is  a  direct 
contrast  to  the  wavelet  which  is  irregular  in  shape.  Its  irregular  shape  improves  the  analysis 
of  signals  with  discontinuities  or  sharp  changes,  while  its  compact  support  allows  temporal 
localization  of  a  signal’s  features  [18]. 

The  two-parameter  family  of  basis  functions  known  as  wavelets  may  be  expressed 

as:  ^a(*)  =  -p  *  [  —  )>  (4-8) 

f  l  a  ) 

where  a  is  a  scale  factor  or  dilation  parameter  and  x  is  a  time  delay.  This  leads  to  the 
definition  of  the  CWT : 

oo  /  \ 

CWTIx.a)  =  —  f  x(»W  —I*.  (4.9) 

fit.  I  4  1 

In  contrast  to  the  Short  Time  Fourier  transform,  the  CWT  offers  improved  frequency 
resolution  at  low  frequency  while  the  time  resolution  is  improved  at  higher  frequency. 
Figure  4.5  illustrates  the  time-frequency  plane  for  the  CWT.  The  wavelet  coefficients, 
determined  by  the  translation  and  dilation  operations  shown  in  Equation  (4.8),  represent  the 
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Repeat  Shifting  Operation 


Figure  4.5:  CWT  Time-Frequency  Display. 


correlation  between  the  wavelet  and  a  localized  section  of  the  signal.  The  process  of 
translation  or  shifting  and  dilation  or  scaling  is  illustrated  in  Figure  4.6. 


Figure  4.6:  Scaling  and  Shifting  Process  of  the  Wavelet  Transform  (Used 
with  permission  from  Joshua  Altmann,  [18]). 
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C.  THE  DISCRETE  WAVELET  TRANSFORM 


The  advantages  of  multi-resolution  analysis  and  the  CWT  have  made  the  wavelet 
transform  an  invaluable  tool  in  signal  processing.  However,  real-time  signal  processing  and 
denoising  operations  using  the  CWT  would  require  a  substantial  expenditure  in 
computational  time  and  resources  [19].  A  much  simpler  implementation  of  the  wavelet 
transform  is  available  which  reduces  the  redundant  information  experienced  in  CWT 
analysis.  Based  on  work  by  Croisier,  Esteban,  and  Galand  [20],  the  discrete  wavelet 
transform  (DWT)  is  an  efficient  discrete  signal  processing  technique  which  lends  itself  to 
digital  computational  methods.  The  DWT  efficiently  performs  signal  decomposition  and 
reconstruction. 

A  complete  mathematical  interpretation  of  wavelets  is  based  on  functional  analysis, 
which  is  well  defined  in  reference  [21].  Basic  concepts  from  both  functional  and  multi¬ 
resolution  analysis,  as  discussed  in  reference  [14],  will  be  used  to  define  the  scaling  function 
and  the  wavelet  function  and  describe  their  significance  in  the  DWT  of  a  signal.  The 
objective,  expressed  by  (4.1),  is  to  represent  a  signal  in  a  given  vector  space  as  a  linear 
combination  of  basis  functions  in  a  transformed  vector  space.  These  basis  functions  include 
dilated  and  shifted  versions  of  the  orthonormal  scaling  and  wavelet  functions.  A  set  of 
scaling  functions  which  span  an  arbitrary  vector  subspace  V0  may  be  defined  as 

(j)*(r)  =  ke  Z,  (j>eL2,  (4.10) 

where  Z  denotes  the  set  of  all  integers  and  L2  is  the  space  of  square  integrable  functions. 
Equation  (4.1)  may  be  rewritten  as, 
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(4.11) 


x(t)  =  YaM)  for  any  x(t)  e  VQ. 

k 

The  subspace  spanned  by  the  scaling  function  may  be  expanded  to  Vj  by  forming  a  set  of 
scaling  functions  which  are  both  scaled  and  translated  as  follows 

=  2^/2<J)(2- h-k).  (4.12) 

Consequently,  (4.11)  becomes 

x(t)  =  Y,  a$j,k^  for  any  x(t)  €  Vf  (4.13) 

k 

For  j  >  0,  (J )jk  is  shorter  and  translates  in  smaller  steps  allowing  the  representation  of  signals 
with  finer  details. 

A  nested  set  of  spanned  spaces  is  defined  as 

-  cV  2cV  jcVqcVjcV^  -  C  L2  (4.14) 

or  Vj  c  Vj+]  for  all  j  e  Z.  (4.15) 

The  recursive  relationship 

4>(t)  =  n  6  z  (4-16) 

n 

allows  the  lower  scale  <j)(t)  to  be  expressed  in  terms  of  a  weighted  sum  of  the  shifted  higher 
scale  (j)(2 1).  The  coefficients  h(n)  are  the  scaling  function  coefficients  which  form  a  filter  to 
be  used  in  the  DWT. 

Considering  the  scaling  function  as  a  coarse  approximation  of  a  signal,  the  wavelet 
function,  which  spans  the  differences  between  the  spaces  spanned  by  the  various  scales  of 
the  scaling  function,  provides  the  finer  details.  The  two  functions  are  orthogonal  at  a  given 
scale  and  the  wavelet  function  can  be  expressed  in  terms  of  the  weighted  sum  of  shifted 
scaling  function  (J)(2 1)  defined  in  (4.16)  by 
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IKO  =  E^V^C2*-")’  n  6  Z,  (4.17) 

n 

where  g(n)  are  the  wavelet  function  coefficients.  The  scaling  function  coefficients  and 
wavelet  function  coefficients  are  required  by  orthogonality  to  be  related  by 

g(n)  =  (-1)"  h(N-l-n),  (4.18) 

where  N  is  the  finite  even  length  of  the  signal.  The  relationship 

V')  =  2JI2¥2jt-k),  (4.19) 

which  is  similar  to  (4.12),  determines  the  two-dimensional  family  of  wavelet  functions  by 
scaling  and  translating  the  mother  wavelet  given  in  (4.17).  Figure  4.7  below  depicts  the 
relationship  between  the  scaling  function  and  wavelet  function  spaces. 

With  the  scaling  and  wavelet  functions  defined,  any  signal  x(t)  €  L2  (R)  can  be  written 


Figure  4.7:  Scaling  function  and  wavelet  function  space 
orthogonal  relationship  [14]. 

as  the  discrete  series 

oo 

x(t)  =  E  c/o(A:)2/o/2(K2 h-k)  +E  E  dfk)2^(2H-k).  (4.20) 

*  *  ;'=;'o 

The  choice  of  j0  sets  the  coarsest  scale  spanned  by  the  scaling  function  providing  a 
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coarse  approximation  of  the  signal.  The  rest  of  the  space  is  spanned  by  the  wavelet  function 
which  provides  the  high  resolution  details  of  the  signal.  The  set  of  coefficients  from  (4.20) 
is  called  the  discrete  wavelet  transform  of  the  signal. 

D.  MULTIRATE  SYSTEMS  ANALYSIS 

Multirate  systems  and  filter  banks  are  the  workhorse  of  the  DWT.  A  deeper 
understanding  of  their  role  in  signal  analysis  and  reconstruction  is  essential  in  exploring  the 
effects  of  combining  the  Wiener  filter  with  wavelet  analysis.  Multirate  system  components 
used  to  implement  the  DWT  include  decimators,  interpolators,  analysis  filters,  and  synthesis 
filters.  These  operations  will  be  discussed  in  the  context  of  the  DWT  and  then  a  simple 
multirate  system  will  be  analyzed. 

Decimation  involves  the  subsampling  or  downsampling  of  a  discrete  sequence  by  a 
factor  of  two  for  the  DWT.  This  is  equivalent  to  discarding  every  other  sample  and  results 
in  reducing  the  sampling  rate  by  a  factor  of  two.  Mathematically  this  process  is  represented 

by 

xl2(n )  =  x(2  n).  (4.21) 

Taking  the  Fourier  transform  of  (4.21)  yields  the  frequency  domain  expression  for 
decimation 

X12(g>)  =  i[X(|)  +  Xq+n)].  (4.22) 

The  Z-transform  of  (4.21)  provides  the  Z-domain  expression 

X,2(z)  -  j[Xfe1>X(-z1'2)].  (4.23) 

The  frequency  and  Z-domain  expressions  show  that  in  addition  to  reproducing  the  input  at 
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half  the  sampling  rate,  X(— ),  decimation  produces  a  second  term  shifted  by  it  radians, 
X(y+Tt).  This  second  term  is  responsible  for  aliasing  which  causes  a  distorted  frequency 
response  unless  the  original  signal  is  properly  bandlimited.  If  X(ca)  is  bandlimited 
to  |  co  |  <  ,  then  there  is  no  overlap  between  the  two  terms;  and  the  aliasing  term  (i.e. 

shifted  term)  can  be  removed  by  filtering. 

Interpolation,  also  known  as  expansion,  upsamples  a  discrete  sequence,  again  by  a 
factor  of  two  for  the  DWT.  The  insertion  of  a  zero  between  each  sample  in  the  sequence 


leads  to  a  doubling  of  the  sampling  rate.  Mathematically  this  is  expressed  as 

’  .  _  f  x(n/2)  for  even  n 

xn\n)  1  o  for  odd  n 


(4.24) 


The  Fourier  and  Z-transforms  of  (4.24)  produce  the  following  results 

X|2(co)  =  X(2u) 

Xn(z)  =  X(z2). 


(4.25) 


Doubling  of  the  sampling  rate  causes  the  original  frequency  response  to  be  compressed  by 


a  factor  of  two.  Additionally,  interpolation  causes  an  effect  known  as  imaging  where  one 
input  frequency  causes  two  output  frequencies,  one  at  radians  and  a  second  image  at 
-j+Tt  radians.  This  is  different  from  aliasing  where  two  input  frequencies,  to  and  ca  +  n  result 
in  the  same  output  [22].  Block  diagrams  of  the  decimation  and  interpolation  operators  are 


Interpolator 


Figure  4.8:  Decimation  and 
Interpolation  Operators. 
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shown  in  Figure  4.8. 

The  analysis/synthesis  filter  bank  employs  the  decimation  and  interpolation  operators 
in  sequential  order  in  order  to  recover  the  original  sequence,  however,  the  effects  of  aliasing 
and  imaging  must  be  considered.  Decimation  followed  by  interpolation  produce  the 
following  frequency  and  Z-domain  results 


X(12)(i2)(g>)  =  1[X(C0)+X(C0+7I)] 

X{lW2)(z)  =  I[X(z)+X(-z)]; 


(4.26) 


which  represent  a  scaled  reproduction  of  the  original  sequence  plus  an  additional  term  caused 
by  aliasing  and  imaging  [22], 

A  lowpass  filter  H0( z)  and  highpass  filter  Hx{ z),  known  as  analysis  filters,  are  now 
introduced  prior  to  decimation  in  order  to  bandlimit  the  input  signal.  These  maximally 
decimated  FIR  filters  form  a  Quadrature  Mirror  Filter  (QMF)  bank  in  which  the  filters 
exhibit  frequency  responses  which  are  a  mirror  image  of  each  other  as  shown  in  Figure  4.9. 
Since  the  analysis  filters  have  a  nonzero  transition  bandwidth  and  stop-band  gain,  some 


Quadrature  Mirror  Filters:  Daubechles— 20  Wavelet 


Figure  4.9:  QMF  bank:  Analysis  Filters  H0  and  Hv 
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aliasing  does  occur  due  to  the  decimation  and  must  be  eliminated  by  other  means. 

The  imaging  problem  is  remedied  by  a  second  QMF  bank  following  the  interpolation 
operation.  These  filters,  known  as  synthesis  filters,  remove  the  images  produced  by  the 
interpolation  operation.  However,  they  suffer  from  the  same  non-ideal  characteristics  as  the 
analysis  filters  and  thus  are  not  capable  of  completely  removing  all  images.  If  the  analysis 
and  synthesis  filters  meet  specific  conditions,  they  are  able  to  remove  all  aliasing  and 
imaging  and  perfectly  reconstruct  the  signal. 

The  analysis  and  synthesis  filters  as  well  as  the  decimator  and  interpolator  operators 
are  now  assembled  to  form  a  filter  bank  as  shown  below  in  Figure  4.10.  Perfect 
reconstruction  is  defined  mathematically  as 

x(n)  =  cx(n-n0 )  (4.27) 


Input  Analysis  Decimators  Interpolators  Synthesis  Output 

Figure  4.10:  Two-channel  Quadrature  Mirror  Filter  bank  [23]. 


where  c  is  a  non-zero  constant  and  n0  a  positive  integer.  Thus  x(n)  is  merely  a  scaled  and 
delayed  version  of  x(n)  [23].  A  Z-transform  analysis  of  the  filter  bank  in  Figure  4.10 
provides  the  conditions  necessary  for  perfect  reconstruction.  The  output  of  the  analysis  filter 


33 


for  channel  k  (k  =  0, 1)  is  expressed  as 


Yt(z)  =  HAz)X(z). 


(4.28) 


According  to  (4.23)  decimation  produces  the  result 

Vk(z)  =  \[Yk(zm)+Yk(-zm)} 

=  ±[Hk(z  m)X(z  m)  +Hk(  -z  m)X(  -z  -m)l 
The  interpolator  output,  determined  by  applying  (4.25)  and  (4.26),  becomes 

Uk(z)  =  Vk(z 2) 

=  i[ HAz)X(z)+Hk(-z)X(-z)]. 


(4.29) 


(4.30) 


X{z)  =  \[X(z)  X(-z)] 


Finally,  the  reconstructed  signal,  which  is  the  sum  of  the  synthesis  filter  bank  channels,  is 

x(Z)  =  ^F0m0(z)*Flmfa))m*^F0m0(-z)*F,m1(-z))x(.-z).  <4.3i) 

In  matrix  form,  (4.31)  becomes 

f»nfe)  »,fe)fF.d 

L0;;.  (4.32) 

-z)  i(z)J 

The  second  component  in  (4.31)  is  due  to  aliasing  and  imaging  and  produces  replica  of  the 
first  term  shifted  ir  radians  to  the  right  on  the  unit  circle.  It  is  commonly  referred  to  as  the 
alias  term  [23].  The  following  must  be  true  to  eliminate  aliasing  and  imaging 

F0(z)tf0(  -z)  +Fl(z)Hl(  -z)  =0.  (4.33) 

Additionally,  the  remaining  term  must  conform  to  the  following  to  remove  distortion  caused 
by  the  presence  of  the  analysis  and  synthesis  filters 

F0(z)Hq(z)+Fx{z)Hx{z)  =  2z _/.  (4.34) 

Having  met  the  alias  cancellation  and  no  distortion  conditions,  perfect  reconstruction  results 
in  the  following,  where  /  is  the  delay  for  this  filter  bank  based  on  the  filter  length: 

X(z)  =  z-‘X(z)  (4.35) 

The  analysis  and  synthesis  filters  can  now  be  selected  based  on  (4.33)  and  (4.34). 
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N 

Assuming  the  following  lowpass  analysis  filter  H0(z)  =  X!  hQ{n)z ,  the  remaining  filters 

n-  0 

are  obtained  as  follows  to  produce  a  two-channel  perfect  reconstruction  filter  bank  [23] 

Hx{z)  =  -z-nH0{-z-x) 

F0(z)  =  z'NH0(Z-1)  (4.36) 

Ffo)  =  z-%(z-1). 

E.  DAUBECHIES  WAVELET 

The  Daubechies  wavelet,  named  after  Ingrid  Daubechies  of  Bell  Laboratories,  is  a 
common  choice  for  the  analysis  and  synthesis  filters,  because  it  possesses  several  nice 
properties.  First,  the  filters  produced  by  this  family  of  wavelets  are  orthogonal  with  compact 
support  providing  many  advantages  previously  mentioned.  Second,  the  frequency  response 
has  maximum  flatness  at  co  =  0  and  co  =  it.  For  a  filter  of  length  N,  placing  half  the  filter 
zeros  at  z  =  - 1  (or  co  =  tc),  results  in  a  maximum  flat  filter  of  regularity  K  =  N/2.  This  second 
property  leads  to  excellent  results  when  Daubechies  filters  are  used  for  the  DWT 
decomposition  and  reconstruction  of  a  large  class  of  signals. 

The  following  theorem  from  Daubechies  1988  paper  on  orthonormal  wavelets  [24] 
summarizes  her  findings: 

Theorem  1.  A  discrete-time  Fourier  transform  of  the  filter  coefficients  h(n )  having 
K  zeros  at  co  =  tc  of  the  form 

H{  co)  =  (lLlilj*C(u)  (4.37) 

satisfies  ltf(co)l2  +  l//(co+7t)l2  =  2  (4.38) 

if  and  only  if  L(co)  =  lC(co)l2  can  be  written  as 
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L(o))  =  P(sin2  o>/2), 

(4.39) 

with 

P(y)  =  PK(y)+yf:R(V2-y), 

(4.40) 

K- 1  ,  > 

where 

W  -  E  T )>* 

\  K  1 

(4.41) 

and  R  is  an  odd  polynomial,  chosen  such  that  P(y)  >0  for  0<y<  1.  There  is  no  explicit 
expression  for  the  scaling  or  wavelet  functions. 

A  brief  explanation  of  Daubechies  development  from  references  [14]  and  [22] 
follows.  Equation  (4.37)  is  expressed  in  terms  of  the  frequency  response  lflr(o))l2  as  follows: 


i#(o>ri  = 


1  +e 


1(0)), 


(4.42) 


where  L(o>)  =  l£(o))l2.  In  trigonometric  terms,  (4.42)  may  be  expressed  as  polynomials  in 
cos(co): 


l//(sin2(co/2))l2  =  [cos2(co/2)F  P(sin2(oo/2))  (4.43) 

which,  after  the  change  of  variables  of  y  =  sin2(co/2)=l-cos2(o)/2),  becomes 

\H(y)\2  =  (1  -y)KP(y),  (4.44) 

where  P(y)  is  an  (N-K)  order  polynomial.  A  similar  derivation  for  I/7(olH-tt)I2  enables  (4.38) 
to  be  expressed  as: 

(l-yfP^+y^l-y)  =  2.  (4.45) 

Daubechies  then  applied  Bezout’s  theorem  with  R= 0  in  (4.40)  resulting  in  minimum  length 
N  for  a  given  regularity  K  (N=K/2)  to  determine  the  explicit  solution  given  by  (4.41).  P(y) 
is  the  binomial  series  for  (1-y)'*,  truncated  after  K  terms,  with  degree  K- 1  and  K  coefficients: 

PK(y)  =  l+Ky+^-y2+-^~fjyK-1  =  (1  -yrK+0(y  K).  (4.46) 


36 


(4.47) 


Hence  the  response  in  terms  of  y  becomes: 

\H(y)\2  =  (i  -y)kY,{K'l+k)yk- 
The  filter  frequency  response  is  now  expressed  as: 

K~\ 

li/(co)l2  =  |1~C08<a>)|t  (4.48) 

where  the  following  transformations  were  used: 

y  =  and  1  -y  =  i22«2>.  (4.49) 

Equation(4.48)  has  K  zeros  at  (o=7t  and  K- 1  zeros  due  to  the  binomial  term  P  (<u). 

Transforming  (4.48)  to  the  Z-domain  using  z  +  zl  =  2  cos(w)  simplifies  the  calculations. 

The  minimum  phase  transfer  function  H(z)  and  its  transpose  H(-z)  are  defined  as  follows: 

N  N 

H(z)  =  £  h(n)z  -n  and  H{z~l)  =  '£h(n)zn.  (4.50) 

n= 0  n=0 

Equation(4.48)  becomes: 

F(Z) = HWHfc")  -  (4-5l) 

which  has  2 K  zeros  at  z  =  -1,  half  of  which  belong  to  H(z),  and  2K-2  zeros  due  to  P(z),  half 
of  which  are  inside  the  unit  circle  and  belong  to  H(z).  The  maximum  flat  Daubechies  filter 
coefficients  are  solved  for  by  applying  spectral  factorization  methods  to  (4.51)  to  obtain  H{z). 
References  [14]  and  [25]  provide  software  methods  for  determining  the  scaling  and  wavelet 
functions  from  which  the  Daubechies  wavelet  filter  coefficients  can  be  determined.  The 
results  using  MATLAB  [25]  are  shown  for  the  Daubechies-2  and  Daubechies-20  wavelet 
(filter  lengths  of  four  and  forty  respectively)  in  Figures  4. 1 1  and  4. 12.  Notice  the  improved 
regularity  and  smoothness  as  the  filter  order  increases. 

The  DWT  is  implemented  through  a  series  of  half-band  highpass  and  lowpass  filters 
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of  varying  cutoff  frequencies  which  analyze  the  signal  at  different  frequency  resolution.  The 
filtering  operations,  which  repeatedly  half  the  signal  spectral  content,  determine  the 
frequency  resolution.  Following  this  filtering  operation,  the  number  of  signal  samples 
remains  unchanged  and  consequently  half  the  samples  are  redundant  based  on  Nyquist’s 
sampling  criterion.  Thus,  the  filtering  operation  is  followed  by  a  downsampling  operation 
to  eliminate  the  unnecessary  samples.  Downsampling  a  signal  by  a  factor  of  two  corresponds 
to  reducing  the  sampling  rate  by  a  factor  of  two  which  results  in  halving  the  time  resolution. 
Thus  each  successive  filtering  and  downsampling  operation  improves  the  frequency 
resolution  by  a  factor  of  two  while  the  time  resolution  suffers  by  a  factor  of  two. 

The  DWT  of  signal  x[n]  is  determined  as  follows.  The  signal  is  passed  through  a 
half  band  digital  lowpass  filter  with  impulse  response  h[n ]  and  a  half  band  digital  highpass 
filter  with  impulse  response  g[n].  The  highpass  filter  is  associated  with  the  wavelet  function 
while  the  lowpass  filter  is  determined  from  the  scaling  function.  This  is  equivalent  to 
convolving  the  signal  with  the  filter  impulse  response  and  is  expressed  as: 

ylow[n ]  =  x[n]*h[n]  =  J^x[k]h[n-k]  =  £ 

Jl  JL  (A  42s) 

yhigh[n]  =  x[n]*g[n]  =  £  4*1  $[»-*]  = 

The  filtering  process  divides  the  frequency  band  of  y[n]  in  half  forming  two  signals,  yhw[n] 
and  yhigl\n\  which  now  have  redundant  information  based  on  Nyquist’s  sampling  criterion. 
These  signals  are  then  decimated  or  subsampled  by  a  factor  of  two  without  any  loss  of 
information.  The  combined  filtering  and  subsampling  process  reduces  the  time  resolution 


38 


Scaling  function  phi 


Wavelet  function  psi 


Figure  4.1 1:  Daubechies-2  Scaling  and  Wavelet  Functions  [25]. 


Scaling  function  phi 


Wavelet  function  psi 
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Analysis  high-pass  filter 


Figure  4.12:  Daubechies-20  Scaling  and  Wavelet  Functions  [25]. 
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by  a  factor  of  two  and  doubles  the  scale  which  improves  the  frequency  resolution  by  a  factor 
of  two.  Mathematically  this  is  expressed  as 

ylow[k]  =  Y,x[n\h[2k-n] 

yhighm  =  1£x[n]g[2k-n].  (4-43) 

n 

The  combined  operation  of  filtering  and  subsampling  is  repeated  to  decompose  the  signal 
to  some  desired  level  which  is  limited  by  the  length  of  the  signal.  At  every  level,  filtering 
and  subsampling  operations  produce  data  with  half  the  number  of  samples  (reducing  the  time 
resolution  by  a  factor  of  two)  and  half  the  frequency  spectrum  bandwidth  (which  doubles  the 
frequency  resolution).  The  filterbank  representation  of  the  DWT  depicted  in  Figure  4.13  can 
be  efficiently  implemented,  and  with  a  reasonable  amount  of  computational  effort  most 
acoustic  signals  can  be  decomposed  for  denoising  purposes. 

The  DWT  coefficients  resulting  from  the  process  described  above  provide  a  spectral 
analysis  of  the  signal  with  a  time  resolution  that  varies  depending  on  the  level  of 
decomposition.  Large  amplitude  DWT  coefficients  indicate  significant  spectral  content  and 
the  position  of  these  coefficients  within  the  DWT  vector  provides  time  localization.  The 
time  resolution  is  precise  at  high  frequencies  and  gets  worse  at  each  successive  level  of 
decomposition  (while  the  frequency  resolution  improves). 
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Level  3  DWT 
Coefficients 


Figure  4.13:  Filterbank  Representation  of  the  DWT  [18]. 


42 


V.  DENOISING  METHODS  USING  WAVELET  THRESHOLDING 


David  Donoho  and  his  colleagues  at  Stanford  University  have  led  the  way  in  applying 
wavelet  techniques  to  the  theory  of  statistics.  A  particularly  relevant  outcome  of  this  work 
is  a  denoising  technique  termed  “wavelet  shrinkage”  by  Donoho  and  his  co-developer  Iain 
Johnstone  [26].  The  wavelet  transform  is  applied  to  a  noisy  data  set  producing  a  set  of 
coefficients  which  are  then  “shrunk”  towards  zero  using  a  soft  or  hard  statistical  thresholding 
method  to  select  the  appropriate  coefficients  which  represent  the  desired  signal.  The 
resulting  coefficients  are  then  back  transformed  to  the  time  domain  to  produce  a  denoised 
signal.  A  brief  and  informal  explanation  of  the  process  follows.  A  more  detailed  argument 
can  be  found  in  many  of  Donoho’ s  writings  [26,27,28]. 

A.  WAVELET  COEFFICIENT  THRESHOLDING 

Assuming  a  suitable  wavelet  basis  function  is  chosen,  the  DWT  decomposition  of 
a  signal  will  compress  the  energy  of  the  signal  into  a  small  number  of  large  magnitude 
wavelet  coefficients.  Gaussian  white  noise  in  any  one  orthogonal  basis  is  transformed  by  the 
DWT  into  white  wavelet  coefficients  of  small  magnitude.  This  property  of  the  DWT  allows 
the  suppression  of  noise  by  applying  a  threshold  which  retains  wavelet  coefficients 
representing  the  signal  and  removes  low  magnitude  coefficients  which  represent  noise. 

Assume  a  finite  signal  s(k)  of  length  N  is  corrupted  by  zero  mean,  additive  white 
Gaussian  noise  n(k)  with  standard  deviation  a,  leading  to  the  noisy  signal 

x(k)  =  s(k)  +  o  n(k).  (5.1) 

The  objective  is  to  recover  the  signal  s(k)  from  the  noisy  observations  x(k)  using  thresholding 


43 


of  the  DWT  coefficients.  The  DWT  of  x(k)  is  performed  using  (4. 18)  and  the  estimate  s  is 

determined  by  thresholding  the  individual  wavelet  coefficients. 

The  three  step  method  for  recovery  of  the  desired  signal,  s(k),  from  the  noisy  data, 
x(k),  follows 

( 1 )  Perform  a  J-level  DWT  of  the  data  yielding  noisy  wavelet  coefficients  Wj  k ,  j=j0, 
...,J ,  k=  0, ... ,  2'  -  1,  where  j  is  the  scale  or  decomposition  level,  and  k  is  the  length 
of  the  signal. 

(2)  Apply  thresholding  in  the  wavelet  domain,  using  either  hard-thresholding  of  the 
coefficients  defined  as: 


™j,k  =  TharlWhk&  = 


Ivv^l  >  6 
lw.  J  <  6  ’ 


(5.2) 


or  soft-thresholding  of  the  coefficients  defined  as: 


siga(wjk)(\wjk\-d) 

0 


\wjJc\  >  6 
1  wjJ  *  8  ‘ 


(5.3) 


The  threshold,  6,  is  usually  determined  in  one  of  four  ways  to  be  described  below. 

(3)  Perform  the  inverse  DWT  and  obtain  the  signal  estimate  s(k). 

Thus,  thresholding  cancels  the  wavelet  coefficients  which  are  below  a  certain 
threshold  value  defined  in  terms  of  the  noise  standard  deviation  o.  The  artificial  signals  used 
in  this  thesis  employ  the  basic  model  (5.1)  with  the  standard  deviation  set  at  one.  The  signal 
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is  then  scaled  to  provide  the  desired  SNR.  In  practice  the  noise  standard  deviation  must  be 
estimated.  Standard  statistical  procedures  estimate  o  as  the  median  absolute  deviation  of  the 
wavelet  coefficients  at  the  finest  level  J,  divided  by  0.6745  [8].  Results  have  shown  that  the 
empirical  wavelet  coefficients  at  the  finest  scale  are,  with  few  exceptions,  essentially  pure 
noise.  Although  subject  to  an  upward  bias  due  to  the  presence  of  some  signal  at  the  finest 
scale,  this  robust  estimate  has  produced  reliable  results  [26].  A  level  dependent 
determination  of  a  is  another  alternative  available  if  colored  noise  is  suspected.  All  three 
options  for  determining  the  noise  standard  deviation  are  available  in  the  MATLAB  Wavelet 
Toolbox  [25], 

B.  THRESHOLD  SELECTION 

Four  commonly  used  threshold  determination  methods  are  discussed  next:  Universal 
threshold,  Stein’s  Unbiased  Risk  Estimator  (SURE),  Hybrid  threshold,  and  the  Minimax 
threshold. 

1.  Universal  Threshold 

The  universal  threshold,  defined  as 

=  oj2fo£N,  (5.4) 

where  a  is  the  standard  deviation  of  the  noise  and  N  is  the  sample  length.  It  uses  a  single 
threshold  for  all  detail  wavelet  coefficients  and  ensures  that  asymptotically  all  detail 
coefficients  are  annihilated  [29].  The  universal  threshold  is  a  less  conservative  method 
which  removes  noise  effectively  but  may  remove  small  details  of  the  signal  which  lie  in  the 
noise  range  [25].  It  is  easily  implemented  and  is  particularly  well  suited  to  minimal  spectral 
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content  signals  which  have  sparse  wavelet  coefficients  [26],  One  drawback  involves  signals 
of  short  duration  (several  hundred  samples  in  length),  in  which  case  other  methods  may 
provide  more  accurate  results  in  a  MSE  sense. 

2.  SURE  Threshold 

The  SURE  threshold,  6S  ,  selects  a  threshold  based  on  Charles  Stein’s  work  on 
unbiased  risk  estimators  [30].  First  implemented  by  Donoho  [27],  this  method  minimizes 
a  risk  function  to  determine  an  optimum  threshold.  The  unbiased  risk  estimate  used  for  soft 
thresholding  is  defined  as  follows: 

RISK(X)  =  1  +  (X2  -  2)  7(1X1  <6,)  +  6*  7(1X1  >  6),  (5.5) 

where  X  ~  7V(0,1),  7(*)  is  the  indicator  function,  and  Ss  is  the  threshold  [29].  Using  the  model 
defined  by  (5.1)  the  threshold  becomes 

6s  =  minj  J2  £  R7Sx(  J,  (5.6) 

where  the  wavelet  coefficients  wjk,  divided  by  the  standard  deviation  have  replaced  the 
variable  X  in  the  unbiased  risk  estimate.  This  technique  suffers  from  the  sparse  wavelet 
coefficient  problem  mentioned  for  the  universal  threshold  and  is  often  combined  with  the 
universal  threshold  in  a  hybrid  threshold  method. 

3.  Hybrid  Threshold 

The  hybrid  threshold  method  is  a  combination  of  the  first  two  methods.  Developed 
by  Donoho  [27],  this  method  uses  the  SURE  threshold  unless  a  low  SNR  situation  with 
sparse  wavelet  coefficients  develops,  in  which  case,  the  universal  threshold  method  is 
utilized. 
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4.  Minimax  Principle  Threshold 

The  minimax  principle,  often  used  in  statistics  to  design  estimators,  uses  a  fixed 
threshold,  6M ,  selected  to  produce  minimax  performance  in  terms  of  the  mean  square  error 
when  compared  to  an  ideal  procedure  [31],  Like  the  universal  threshold,  minimax  uses  a 
single  threshold  for  all  detail  wavelet  coefficients.  The  minimax  threshold  is  a  function  of 
the  sample  size,  the  threshold  type  (soft  or  hard),  and  the  oracle  type  (projection  or 
shrinkage)  and  can  be  found  in  tabulated  form  in  references  [28]  and  [29].  MATLAB 
Wavelet  Toolbox  [25]  uses  the  following  expression  which  approximates  the  minimax 
threshold  for  the  soft  threshold  nonlinearity,  with  comparison  to  a  projection  oracle: 

6  Jn)  =  0.3936  +  0.1829*(log(«)/log(2)),  (5.7) 

where  n  is  the  sample  length. 

5.  Threshold  Discussion 

The  four  thresholding  methods  are  compared  in  Figure  5.1  where  various  wavelet 
thresholding  schemes  are  applied  to  the  same  noisy  sinusoidal  signal  used  earlier  in  Chapter 
HI.  The  SURE  and  Hybrid  thresholding  methods  represent  the  signal  as  a  Smoothed  sinusoid 
after  a  four-level  decomposition  using  the  Daubechies-20  wavelet  for  the  DWT.  The 
Universal  and  Minimax  thresholding  methods  applied  to  a  four-level  decomposition  do  not 
perform  as  well  as  expected,  although  this  signal  has  minimal  spectral  content  and  therefore 
a  sparse  representation  in  the  wavelet  domain  as  seen  in  Figure  5.2.  This  is  explained  by 
comparing  Figures  5.2, 5.3,  and  5.4  which  show  the  approximation  and  detail  coefficients 
before  and  after  universal  and  Minimax  thresholding.  The  Universal  and  Minimax  methods 
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apply  a  constant  threshold  for  a  given  decomposition  level.  This  results  in  poor  denoising 
performance  at  any  scale  other  than  the  coarsest  or  highest  scale  otherwise  known  as  the 
approximation.  This  is  true  because  the  high  threshold  value  results  in  the  signal  as  well  as 
the  noise  being  removed  by  the  thresholding  except  for  the  approximation  which  is  not 
thresholded  in  “wave  shrink”  methods.  This  effect  is  most  noticeable  at  low  SNR  values 
where  the  signal  and  noise  wavelet  coefficients  approach  the  same  magnitude.  The  SURE 
and  Hybrid  thresholding  methods  employ  scale  dependent  thresholds  which  are  more 
effective  in  a  low  SNR  environment  as  demonstrated  in  Figures  5.5  and  5.6  which  show  the 
wavelet  coefficients  after  applying  these  two  methods. 

The  normalized  spectral  decomposition  filters  are  pictured  in  Figure  5.7  for  the 
Daubechies-20  wavelet.  A  MSE  comparison  of  the  four  thresholding  methods  across  the 
normalized  frequency  spectrum  from  0  to  0.5  (Sampling  Frequency  =  1)  is  shown  in  Figures 
5.8  and  5.9.  The  poor  performance  of  the  Universal  and  Minimax  methods  at  low  SNR  level 
is  clearly  evident.  Another  significant  phenomenon  noticeable  in  these  figures  is  the 
improvement  in  normalized  MSE  with  each  successive  scale  and  decomposition  level.  This 
occurs  because  with  each  successive  scale  there  is  less  noise  energy  as  represented  by  the 
wavelet  coefficients,  whereas  the  signal  energy  remains  the  same.  The  spectral  bandwidth 
is  reduced  with  each  level  of  decomposition  eliminating  approximately  half  the  remaining 
white  noise  which  reduces  the  noise  variance  and  improves  the  SNR  at  that  particular  scale. 
The  energy  difference  between  the  signal  and  noise  wavelet  coefficients  is  greater  at  higher 
or  coarser  scales  resulting  in  more  effective  removal  of  the  noise  by  thresholding. 
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A  closer  inspection  of  the  threshold  values  chosen  for  a  particular  thresholding 
method  reveals  that  the  threshold  is  dependent  on  the  number  of  decomposition  levels  even 
though  the  signal  length  of  1024  data  point  remains  the  same.  This  behavior  is  due  to  the 
specific  implementation  employed  in  the  MATLAB  Wavelet  Toolbox  [25]  which  doesn’t 
keep  the  total  number  of  wavelet  coefficients  constant.  Threshold  values  for  the  Universal 
and  Minimax  methods  are  compared  below  for  four  levels  of  decomposition. 


THRESHOLD  VALUES  ( SIGNAL  LENGTH  =  1024  SAMPLES) 

Decomposition 

Level 

Total  Number  of 
Wavelet  Coefficients 

Universal 

Threshold 

Minimax 

Threshold 

1 

1062 

3.7331 

2.2322 

2 

1101 

3.7427 

2.2417 

3 

1140 

3.7520 

2.2509 

4 

1178 

3.7607 

2.2596 

Since  the  threshold  value  for  the  Universal  and  Minimax  methods  is  dependent  on  the 
number  of  wavelet  coefficients,  the  threshold  varies  depending  on  the  decomposition  level 
resulting  in  independent  MSE  curves.  The  SURE  and  Hybrid  methods  calculate  a  different 
threshold  for  each  scale  allowing  them  to  generally  outperform  the  Universal  and  Minimax 
methods.  This  is  particularly  noticeable  at  the  lower  or  finest  scale  where  scale  dependent 
SURE  and  Hybrid  thresholding  methods  generate  much  better  results. 
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Original  Signal  and  Noisy  Signal:  SNR  =  OdB,  Signal  Freq  =  .1 


Figure  5.1:  Comparison  of  wavelet  thresholding  schemes;  sinusoidal  signal  at  Frequency 
=0.1,  Sampling  Frequency  Fs=l,  Original  signal  (dashed);  Denoised  signal  (solid);  Five- 
level  DWT  decomposition  using  Daubechies-20  wavelet. 
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Approximation  Coef.:  Freq=0.1,  SNR=OdB 


Figure  5.2:  Approximation  and  detail  wavelet  coefficients  resulting  from  five-level  DWT 
decomposition  of  noisy  sinusoidal  signal,  Frequency=0.1  (Fs=l),  SNR=0  dB  using 
Daubechies-20  wavelet. 


Universal  Threshold,  Approx.  Coef.:Freq=0.1,  SNR=0dB 


Figure  5.3:  Approximation  and  detail  wavelet  coefficients  after  applying  Universal  soft 
thresholding.  Threshold  values  indicated  for  each  level  or  scale. 
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0  100  200  300  400  500  600 

Figure  5.4:  Approximation  and  detail  wavelet  coefficients  after  applying  Minimax  soft 
thresholding.  Threshold  values  indicated  for  each  level  or  scale. 
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Figure  5.5:  Approximation  and  detail  wavelet  coefficients  after  applying  SURE  soft 
thresholding .  Threshold  values  indicated  for  each  level  or  scale. 
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Hybrid  Threshold,  Approx.  Coef.:Freq=0.1,  SNR=0dB 


Figure  5.6:  Approximation  and  detail  wavelet  coefficients  after  applying  Hybrid  soft 
thresholding.  Threshold  values  indicated  at  each  level  or  scale. 
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Figure  5.8:  Soft  thresholding  methods  applied  to  noisy  sinusoidal  signal:  SNR  =  10  dB 
Four  level  DWT  decomposition  using  Daubechies-20  wavelet.  MSE  curves  for  one 
through  four  levels  (Lvl)  of  decomposition  or  scale  pictured.  Ten  trial  average. 


VL  FIR  WIENER  WAVELET  FILTER 


Wiener  filtering  and  wavelet  thresholding  are  both  effective  means  of  denoising 
acoustic  signals.  In  this  chapter  the  two  methods  are  combined  to  form  an  enhanced  FIR 
filter  which  will  be  referred  to  as  the  Wiener  Wavelet  filter.  This  method  of  denoising 
utilizes  the  discrete  wavelet  transform  to  represent  the  signal  as  a  series  of  wavelet  expansion 
coefficients.  Wavelet  analysis  has  several  distinct  properties  which  make  it  an  ideal  method 
of  signal  decomposition,  analysis,  and  reconstruction.  First,  wavelets  are  an  unconditional 
basis  for  a  wide  variety  of  signals  which  means  that  a  signal  can  be  represented  by  a  small 
number  of  expansion  coefficients,  dj  k  in  (4.20),  since  the  magnitudes  of  these  coefficients 
drop  off  rapidly  for  a  wide  class  of  signals  [32].  Second,  wavelet  analysis  provides  a  multi¬ 
resolution  analysis  in  time  and  frequency  allowing  a  more  accurate  description  and 
separation  of  signal  and  noise  [14].  Third,  there  are  a  variety  of  wavelet  functions  which 
makes  wavelets  adaptable  to  represent  signals  of  differing  characteristics.  Finally,  wavelet 
analysis  and  particularly  the  DWT  is  well  suited  to  implementation  on  a  digital  computer 
since  it  involves  only  multiplications  and  additions  and  not  derivatives  or  integrals. 

Once  the  signal  has  been  transformed  into  the  wavelet  domain,  the  wavelet  expansion 
coefficients  are  filtered  using  the  optimal  Wiener  filter.  In  this  process  the  wavelet 
coefficients  are  being  shrunk  towards  zero,  not  by  a  nonlinear  thresholding  method  but  by 
an  optimal  linear  filtering  process.  The  shrunk  coefficients  are  then  inverse  transformed  to 
reconstruct  the  denoised  signal.  The  Wiener  filtering  process  introduces  distortion  and 
aliasing  which  affects  the  ability  to  perfectly  reconstruct  the  denoised  signal.  These  effects 
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will  be  discussed  after  introducing  the  necessary  multirate  system  and  filter  bank  theory. 

A.  MULTIRATE  SYSTEMS  ANALYSIS 

The  application  of  multirate  systems  and  filter  banks  in  signal  analysis  and 
reconstruction  was  discussed  in  Chapter  IV.  The  necessary  conditions  for  perfect 
reconstruction  were  derived.  The  optimal  Wiener  filter  will  now  be  applied  following 
decimation  and  prior  to  interpolation  in  each  stage  of  a  filter  bank.  The  two-channel  filter 
bank  is  illustrated  below  in  Figure  6.1. 


Input  Analysis  Decimators  Wiener  Filters  Interpolators  Synthesis  Output 
Figure  6.1:  FIR  W iener  W avelet  Filter  bank. 


The  Wiener  filters  which  have  been  inserted  in  the  filter  bank  above  must  be  provided 
with  the  statistical  properties  of  the  noise.  This  is  accomplished  by  isolating  a  noise-only 
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portion  of  the  noisy  signal,  x(n).  A  DWT  of  the  noise,  w(n),  is  performed  using  the  same 
analysis  filter  bank  structure,  as  was  used  for  the  noisy  signal.  The  Wiener-Hopf  equations 
(3.6)  to  (3.8)  are  applied  to  the  noisy  signal  wavelet  coefficients  and  the  noise  wavelet 
coefficients  to  produce  the  filter  coefficients  for  the  Wiener  filters  W0  and  Wj.  This  two- 
channel  filter  bank  is  repeated  to  produce  a  multi-level  FIR  Wiener  wavelet  filter  bank. 

A  four-level  decomposition  FIR  Wiener  wavelet  filter  bank  is  applied  to  the  noisy 
sinusoidal  signal  from  previous  chapters  and  the  results  are  displayed  in  Figure  6.2.  The 
normalized  analysis  filters  for  the  Daubechies-20  wavelet  used  in  this  analysis  are  pictured 
in  Figure  6.3.  Figure  6.4  and  6.5  show  the  MSE  performance  for  this  filter  bank  with  the 
Wiener  filter  order  varying  from  4  to  20.  In  Figure  6.4  an  independent  white  Gaussian  noise 
sample  with  a  zero  mean  was  provided  to  the  Wiener  filter  while  in  Figure  6.5  the  actual 
noise  sample  used  to  produce  the  noisy  signal  was  provided.  The  MSE  difference  between 
the  independent  and  actual  noise  cases  becomes  greater  at  lower  SNR  levels.  A  common 
feature  of  all  wavelet  based  techniques  is  the  improved  MSE  performance  with  successive 
levels  of  decomposition.  This  was  observed  in  wavelet  thresholding  and  will  occur  in  the 
HR  Wiener  wavelet  methods  described  in  the  next  chapter.  The  Wiener  filters  are  applied  on 
a  level  dependent  basis,  which  is  different  from  universal  and  minimax  wavelet  thresholding 
methods  where  the  threshold  was  applied  in  a  global  manner  across  all  of  the  scales  of  a 
particular  level  of  decomposition  resulting  in  independent  MSE  curves. 

Another  characteristic  observed  in  wavelet  based  denoising  methods  is  the  degraded 
MSE  performance  in  the  transition  regions  of  the  analysis  and  synthesis  filters.  The  DWT 
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filter  bank,  as  originally  constructed,  met  the  perfect  reconstruction  conditions. 
Consequently,  aliasing  and  distortion  effects  which  would  normally  occur  in  the  transition 
regions  of  the  filters  are  canceled.  Denoising  by  filtering  or  thresholding  the  DWT 
coefficients  between  the  decimation  and  interpolation  operations  results  in  a  filter  bank 
which  no  longer  meets  the  perfect  reconstruction  conditions.  This  effect  is  most  noticeable 
in  the  filter  transition  regions  because  this  is  where  the  perfect  reconstruction  property  was 
so  essential  in  canceling  the  effects  of  aliasing  and  distortion.  This  unfortunate  result  can  be 
minimized  by  choosing  wavelet  analysis  and  synthesis  filters  of  longer  length  to  reduce  the 
transition  region.  Some  detrimental  effects  occur  as  a  result,  including  more  computational 
effort  and  longer  filter  transient  periods.  The  transient  period,  during  which  the  filter  output 
is  unreliable,  is  equal  to  the  filter  length.  As  the  filter  length  increases  relative  to  the  signal 
length  the  transient  period  may  become  significant. 
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Orlo l rtetl  (detshod)  Noisy  Signal:  SNR  —  OcJB.  Signal  Fraq  -■  0.1 


Original  «.r»d  f=IF*  W»  «n«r  Wavalot  Filtarad  Signal:  MS  EE—  0.06392 


Figure  6.2:  FIR  Wiener  Wavelet  Filter  (a)  Original  (dashed)  and  Noisy  Sinusoidal 
Signal.  Frequency=0.1,  Sampling  Frequency=l.  (B)  Original  (dashed)  and  filtered 
signal.  Three  level  decomposition  using  Daubechies-20  wavelet.  Filter  Order  =  8. 


Four  Level  Spectral  Decomposition  Using  Daubechies — 20  Wavelet 


Figure  6.3:  Gain  normalized  analysis  filters  using  Daubechies-20  wavelet. 
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FIR  Wiener  Wavelet:  Order  4 


FIR  Wiener  Wavelet:  Order  8 


FIR  Wiener  Wavelet:  Order  12  FIR  Wiener  Wavelet:  Order  20 


Figure  6.5:  FIR  Wiener  Wavelet  filter  applied  to  noisy  sinusoidal  signal:  SNR  =  0  dB. 
Four  level  DWT  decomposition  using  Daubechies-20  wavelet.  Actual  noise  sample.  Ten 
trial  average. 
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B.  ALIAS  CANCELLATION 

The  conditions  for  perfect  reconstruction  were  specified  in  Chapter  4.  A  two-channel 
filter  bank  was  analyzed  in  the  Z-domain  with  the  following  result 

t(,z)  -  (6.1) 

Inserting  the  Wiener  filters  to  denoise  the  data  in  the  wavelet  domain  alters  (6.1)  as  follows: 

m  = 

i  /  \  (6.2) 

+^F0(z)W0(z2)H0(-z^Fl(z)Wl(.z' 2)Hx(-z))X(-z). 

Alias  cancellation  is  obtained  by  choosing  the  following  synthesis  filters, 

F0(z)  *  W,(z2)H,(-z) 

F,(z)  =  -W0(z2)H0(-z). 

which  cancel  the  X(-z)  alias  term  in  (6.2)  leading  to 

X(z)  =  ±[Wx(z  2)Wq(z  2){Hq{z)Hx{-z)  ~H0(  -z)Hx(z)]]X(z)  (6.4) 

As  before,  the  alternating  flip  relationship  between  the  analysis  filters  is  chosen  to  eliminate 
distortion  affects: 

Hx(z)  =  -z~nH0(-Z-').  (6.5) 

Substituting  (6.5)  into  (6.4)  results  in  the  following  Z-domain  expression  for  the 
reconstructed  signal  in  terms  of  the  original  signal, 

X(z)  =  Wx(z  2)W0(z  2)z  -NX(z).  (6.6) 

Equation  (6.6)  is  applied  to  the  two-level  decomposition  pictured  in  Figure  6.6  with  the 
following  result  for  the  reconstructed  Z-domain  signal, 

X(z)  =  WqCz  2)  Wj (z 4) W2(z  4)z  ~NX(z).  (6.7) 
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(6.8) 


Equations  (6.6)  and  (6.7)  can  be  generalized  for  a  J  level  decomposition  as  follows, 

m  =W0(Z2)Z~NX(Z)  flWk(z2J). 

k= 1 

The  reconstructed  signal  is  now  a  function  of  a  delay  determined  by  filter  length  N,  and  the 
product  of  a  series  of  alias  cancellation  filters  given  in  (6.8).  These  alias  cancellation  filters 
are  nothing  more  than  interpolated  Wiener  Filters  which  are  multiplied  in  the  Z-domain  or 
convolved  in  the  time  domain.  Unfortunately,  simulations  show  that  the  alias  cancellation 
filters  produce  significant  phase  and  amplitude  distortion  which  is  more  detrimental  than  the 


Analysis  Filters  Wiener  Filters  Synthesis  Filters  Alias  Cancellation  Filters 

Figure  6.6:  Two-level  FIR  Wiener  Wavelet  Filter  Bank  with  Alias  Cancellation 

aliasing  effects  they  are  designed  to  remove. 

The  FIR  Wiener  Wavelet  filter  with  alias  cancellation  is  applied  to  the  same  noisy 
sinusoidal  signal  as  used  before.  The  denoised  signal  is  poorly  reconstructed  based  on  an 
evaluation  of  the  MSE  alone.  The  MSE  values  are  high  with  erratic  results  across  the 
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spectrum.  A  comparison  of  the  power  spectral  densities  of  the  noisy  signal  and  the  denoised 
signal  for  two-level  FIR  Wiener  wavelet  filters  with  and  without  alias  cancellation  is  shown 
in  Figure  6.7.  Figures  6.7  (a)  and  (b)  show  the  power  spectral  densities  of  the  noisy  and 
denoised  signal  without  alias  cancellation.  The  alias  cancellation  filter  with  frequency 
response  as  shown  in  Figure  6.7  (c)  is  used  to  produce  the  denoised  signal  shown  in  Figure 
6.7  (d)  which  has  been  normalized  to  remove  the  effects  of  the  gain  attenuation.  Although 
the  aliasing  frequency  at  a  normalized  frequency  of  0.2  has  been  cancelled,  significant 
distortion  is  occuring  at  this  same  frequency  as  well  as  at  normalized  frequencies  of  0.05  and 
0.45.  The  average  gain  attenuation  for  a  two-level  FIR  Wiener  wavelet  decomposition  with 
alias  cancellation  (but  without  energy  normalization  of  signal  spectra)  is  41.5  dB  measured 
at  the  signal  frequency  which  is  quite  significant  when  compared  to  0.49  dB  for  the  same 
filter  without  alias  cancellation  (Attenuation  values  for  signal  frequencies  from  0.01  to  0.49 
were  computed  and  then  averaged). 

The  FIR  Wiener  filter  with  alias  cancellation  might  still  be  considered  a  suitable 
method  if  the  denoised  signal  can  be  normalized.  A  final  comparison  is  made  between  the 
two  filtering  methods  based  on  the  difference  of  two  spectra,  determined  from  energy 
normalized  signals,  on  a  log  magnitude  versus  frequency  scale  defined  by: 

V(o)  =  log  5(o>)  -  log  5/(co).  (6.9) 

where  5(co)is  the  energy  normalized  noise-free  signal  spectra  and  5/(<o)is  the  energy 
normalized  denoised  signal  spectra.  The  Mean  Absolute  Log  Spectrum  Distortion  (MALSD) 
is  defined  in  reference  [33]  as: 
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MALSD  =  |  IV(o>)l^.  (6.10) 

-71 

Expressed  in  terms  of  a  summation  of  discrete  terms,  (6.10)  becomes 

AM 

MALSD  =  Y,  W(ni/N)\.  (6.11) 

i=0 

where  N  is  the  number  of  frequency  samples  taken  between  zero  and  the  half  sampling 
frequency,  7t.  When  applied  to  two-level  decomposition  filter  banks  with  and  without  alias 


cancellation,  the  following  results  were  obtained  (ten  trial  average): 


FIR  Wiener  Wavelet  Filter 

MALSD 

Without  Alias  Cancellation 

1703 

With  Alias  Cancellation 

2584 

The  FIR  Wiener  Wavelet  filter  without  alias  cancellation  exhibits  a  smaller  magnitude 
MALSD  indicating  the  addition  of  alias  cancellation  filters  is  detrimental  to  the  denoising 
performance  of  this  filter.  For  this  reason  any  further  results  obtained  using  a  FIR  Wiener 
Wavelet  filter  will  not  include  alias  cancellation  filters. 
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Amplitude  (dB) 


30 


(a)  Noisy  Signal:  SNR  =  OdB 


(b)  FIR  Wien.  Wave.  No  Alias  Cane. 


Figure  6.7:  Power  spectral  densities  of:  (a)  noisy  sinusoidal  signal,  Frequency=0.3,  Fs=l, 
SNR=0  dB,  (b)  denoised  result  as  determined  with  no  alias  cancellation  filter  applied,  (d) 
denoised  result  with  alias  cancellation  filter  in  (c)  applied. 


VII.  HR  WIENER  WAVELET  FILTER 


The  denoising  methods  developed  thus  far  have  been  restricted  to  the  FIR  Wiener 
filter,  Wavelet  thresholding  methods,  and  a  combination  of  the  two.  This  chapter  develops 
the  noncausal  HR  Wiener  filter  and  then  introduces  its  application  in  the  wavelet  domain. 
Next,  the  results  are  compared  to  previously  discussed  methods. 

A.  HR  WIENER  FILTER 

The  HR  Wiener  filter  is  recursive  in  form  and  requires  fewer  parameters  to  determine 
the  optimal  filter  weights  than  a  comparable  FIR  form  of  the  filter  [12] .  The  HR  causal  filter 
problem  was  solved  (originally  in  the  continuous  case)  by  Wiener  in  the  transform  domain 
using  spectral  factorization  methods  [11].  This  section  describes  the  transform  domain 
solution  of  the  noncausal  Wiener-Hopf  equation.  The  noncausal  solution  allows  the  filter 
to  “look  ahead”  of  real  time  and  thus  operates  in  an  off-line  application.  This  advantage 
results  in  improved  performance  in  the  mean-square  error  sense  over  the  causal  FIR  Wiener 
filter. 

The  noncausal  HR  Wiener  filter  estimate  of  the  desired  signal  is  of  the  form 

00 

s(n)  =  y(n)  =  Y  h  *(k)x(n - k ),  (7.1) 

which  differs  from  (3.1)  in  two  respects.  First,  the  upper  limit  of  the  sum  extends  to  +°°  and 
the  impulse  response  has  nonzero  values  for  n  <  0.  Application  of  the  orthogonality  principle 

yields  the  noncausal  HR  form  of  the  Wiener-Hopf  equation 

00 

£  RJ!-k)h  *(k)  =  R*x(i);  <  i  <  ~  (7.2) 

k=-°° 
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which  can  be  written  in  convolutional  form  as: 

Rx(i)*h\i)  =  Rs'x(i),  (7.3) 

with  mean-square  error 

oo 

o\  =  Rs{Q>)~Y  h*msx(k).  (7.4) 

After  conjugating  equation  (7.3),  the  left-hand  side  is  a  discrete  convolution  of  the  filter 
coefficients  h(k),  and  the  auto-correlation  function,  Rx(k).  Expressing  the  Wiener-Hopf 
equation  in  the  frequency  domain  leads  to: 

Sx(u)H(  to)  =  (7.5) 


where 


H((x))  =  h(n)e  "1C0n,  -7t<co<7t 

n~-°° 

00 

ssx( CO)  =  £  /?  -7t^O)<7t 

n  =  -oo 
oo 

5x(co)  =  Rx{ri)e~iu>n,  _7t<o)<TC. 


(7.6) 


Thus  the  filter  transfer  function  H(oi)  can  be  expressed  as: 


HJo)  « 


SXu) 


(7.7) 


Assuming  that  the  signal  and  noise  are  independent  and  have  a  zero  mean  results  in: 

SJ(o)=S((o)  and  Sj (to)  =St(co)  +Sj (co) ,  (7.8) 


allowing  (7.7)  to  be  expressed  as 

tf„,(0>)  = 


5»  -  S» 

5r(co)  ’ 


(7.9) 
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which  may  be  rewritten  as 


SM/SJp) 

H  .(co)  =  - s~— — —  . 

"c  S  ((ji)/S  (co)  +  1 


(7.10) 


The  expression  Ss(<ui)/Sw(o))  is  a  measure  of  the  signal-to-noise  power  ratio  at 
frequency  co.  When  5i(co)/5M,(co) »  1  the  filter  coefficient  magnitude,  //nc(co),  approaches 


unity,  and  when  Ss((j))/Sw(u>)  « 1  the  filter  coefficient  magnitude  approaches  zero.  Between 
these  two  limits  the  filter  coefficients  are  chosen  to  pass  the  signal  and  eliminate  noise  with 
minimal  distortion  [34]. 

B.  IIR  WIENER  FILTER  APPLIED  TO  THE  WAVELET  DOMAIN 


A  key  property  which  underlies  the  success  of  wavelets  is  that  they  form  an 
unconditional  basis  for  a  wide  variety  of  signal  classes  [35].  Wavelet  expansions  can 
therefore  concentrate  the  signal  energy  into  a  relatively  small  number  of  large  coefficients. 
This  signal  compaction  property  makes  the  discrete  wavelet  transform  an  ideal  tool  in 
constructing  an  empirical  HR  Wiener  filter.  The  result  will  be  referred  to  as  the  HR  Wiener 
Wavelet  filter. 

The  original  time-series  model  given  by 

x(ri)  =  s(ri)  +  w(n);  n  =  1,2  ,.../V 
becomes,  in  the  wavelet  domain, 

yjk  -  0jk  +  Zjk  7=1,2,...,/  ,  k=l,2,,..JV, 

where  the  terms  are  the  DWT  coefficients  found  from  the  following 

x(n)  =  J^yjk2ja^(2jt-k) 
jjc 

s(n)  =  '£Qjk2/2W2jt-k) 

J* 

w(n)  =  2j/Z  ty(2jt~k). 

J* 


(7.11) 

(7.12) 


(7.13) 
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The  indices  j  and  k  refer  to  level  of  DWT  decomposition  and  signal  sample  number 
respectively. 

The  signal  estimate  s  is  determined  by  applying  the  inverse  wavelet  transform  to  the 
coefficients,  0.^.  The  noise  used  to  produce  the  noisy  signal  is  provided  in  order  to 
determine  the  optimal  Wiener  filter  coefficients.  Denoising  an  actual  ocean  acoustic  signal 
would  require  an  isolated  noise  only  segment  or  some  other  noise  estimate  in  order  to 
determine  the  approprite  noise  statistical  characteristics.  A  noise  sample  or  estimate  which 
is  not  highly  correlated  with  the  actual  noise  will  degrade  the  denoising  capability  of  this 
filter.  The  observation  sequence  wavelet  coefficients,  yjk ,  and  noise  wavelet  coefficients, 
zjrk ,  are  provided  allowing  (7.9)  to  be  rewritten  as  follows  for  a  multiple  level  DWT: 


m  m 


m 


N_ 

2j’ 


j  =  1,2,..., 7 


(7.14) 


where  N  is  the  signal  length,  j  is  the  level  of  the  DWT  decomposition,  and  J  represents  the 
maximum  decomposition  level.  This  filter  is  applied  to  the  observation  sequence  wavelet 
coefficients  at  each  decomposition  level  to  determine  the  signal  estimate  coefficients  as 
follows 


j* 


=  hjyjjc,  7=1,2,...,/, 


k= 1,2,.../V2‘/ 


(7.15) 
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The  DR  Wiener  Wavelet  method  is  implemented  in  a  one-level  filter  bank  structure, 
as  illustrated  in  Figure  7.1.  This  is  the  same  filter  bank  structure  as  that  used  for  wavelet 
thresholding  and  FIR  Wiener  wavelet  filtering  except  that  now  at  each  scale  the  noisy 
wavelet  coefficients  are  thresholded  by  the  filter  coefficients  hx  and  .  The  structure  may 
be  expanded  to  additional  levels  or  scales,  limited  only  by  the  signal  length.  This  denoising 
method  is  applied  to  the  noisy  sinusoid  used  in  previous  methods  with  results  as  shown  in 
Figure  7.2.  The  filter  represents  the  denoised  signal  as  a  smoothed  sinusoid  after  a  three 
level  decomposition  using  the  Daubechies-20  wavelet.  The  Daubechies-20  analysis  filters 


Input  Analysis  Decimators  IIR  Wiener  Filters  Interpolators  Synthesis  Output 
Figure  7.1:  One-level  noncausal  HR  Wiener  wavelet  filter  bank, 
of  length  forty  are  pictured  in  Figure  7.3  with  gain  normalized  to  the  maximum  gain  filter 
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in  order  to  show  the  spectral  bands  associated  with  the  scales.  This  normalization  is  only 
used  to  allow  a  better  spectral  view  of  the  filters  and  is  not  representative  of  the  actual  filter 
gains  used  in  the  DWT. 


Original  (dashed)  and  Noisy  Signal:  SNR  =  OdB,  Signal  Freq  =  0.1 


Original  and  HR  Wiener  Wavelet  Filtered  Signal:  MSE  =0.1647 


O  5  10  15  20  25  30  35  40  45  50 


Time  (Sample  Number) 

Figure  7.2:  HR  Wiener  Wavelet  Filter  (a)  Original  (dashed)  and  Noisy 
Sinusoidal  Signal  Frequency=0.1,  Fs=l.  (b)  Original  (dashed)  and 
filtered  signal.  Three  level  decomposition  using  Daubechies-20  wavelet. 

A  ten  trial  average  normalized  MSE  is  plotted  across  the  normalized  frequency  range 
of  0  to  0.5  for  a  sampling  frequency  of  1  as  shown  in  Figures  7.4  and  7.5.  An  independent 
noise  sample  of  white  Gaussian  noise  with  a  zero  mean  was  used  in  Figure  7.4  while  the 
actual  noise  sample  was  used  in  Figure  7.5.  The  MSE  is  reduced  with  each  successive  scale 
or  decomposition  level.  This  occurs  because  with  each  successive  scale  there  is  less  noise 
energy  as  represented  by  the  wavelet  coefficients,  whereas  the  signal  energy  remains  the 
same.  The  spectral  bandwidth  is  reduced  with  each  level  of  decomposition  eliminating 
approximately  half  the  remaining  white  noise  which  reduces  the  noise  variance  and  improves 


the  SNR  at  that  particular  scale.  This  leads  to  an  DR  Wiener  filter  which  is  more  effective 
at  removing  noise  as  demonstrated  by  the  improved  normalized  MSE. 

The  DR  Wiener  filter  is  a  level  dependent  filter  which  means  that  increasing  the  level 
of  decomposition  will  only  improve  the  MSE  performance  for  the  spectral  band  associated 
with  the  particular  analysis  filter.  This  concept  is  better  understood  by  comparing  Figures 
7.3  and  7.4.  Notice  the  degraded  MSE  that  occurs  in  the  analysis  filter  overlap  regions.  The 
spectral  range  over  which  this  degradation  occurs  can  be  reduced  by  choosing  analysis  filters 
with  smaller  transition  regions.  Higher  order  filters  will  accomplish  this  at  the  expense  of 
greater  computational  effort.  The  distortion  can  never  be  completely  removed  since  the 
perfect  reconstruction  property  of  the  filter  bank  has  not  been  met  due  to  the  insertion  of  the 
DR  Wiener  filter. 


Four  Level  Spectral  Decomposition  Using  Daubechies-20  Wavelet 


Figure  7.3:  Normalized  Spectral  Decomposition  Filters  Used  For  DR 
Wiener  Wavelet  Denoising.  (Daubechies-20  Wavelet). 
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Figure  7.5:  Mean-square  error  comparison  for  four  level  decomposition  of  noisy  sinusoid 
using  Daubechies-20  wavelet.  Actual  noise  sample.  Average  of  10  trials. 
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C.  WIENERSHRINK 


The  HR  Wiener  wavelet  filter  just  discussed  requires  an  isolated  noise  segment  in 
order  to  determine  the  statistical  noise  characteristics  necessary  to  generate  a  Wiener  filter. 
WienerShrink,  developed  by  Ghael,  Sayeed,  and  Baraniuk  of  Rice  University  [36],  is  a 
wavelet-based  Wiener  filtering  technique  which  differs  from  previous  methods  in  this  regard. 
Initially,  a  wavelet  thresholding  technique  similar  to  those  discussed  in  Chapter  V  is  applied 
to  the  noisy  signal.  The  resulting  denoised  signal,  referred  to  as  the  pilot  signal,  is  then 
wavelet  transformed  a  second  time  and  used  to  construct  an  HR  Wiener  filter.  This  wavelet- 
based  empirical  Wiener  filtering  method  is  described  by  the  block  diagram  shown  in  Figure 
7.6. 

The  noisy  signal  as  represented  by  (7.11)  is  wavelet  transformed  by  the  DWT 


x  =  s  +  w 

- > 


Figure  7.6:  Wavelet-based  empirical  Wiener  filtering  [25]. 
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represented  by  Wj.  Any  of  the  four  wavelet  thresholding  methods,  represented  by  block  H, 
is  applied  to  the  wavelet  coefficients  ya>^,  to  produce  the  wavelet  coefficient  signal  estimate, 
0(1).  .  This  is  then  transformed  by  the  inverse  DWT  represented  by  Wf1  to  form  the  signal 
estimate,  sx .  The  signal  estimate  is  then  wavelet  transformed  again  in  block  W2  to  form  the 
wavelet  coefficient  signal  estimate  0(21)jr  This  is  used  to  design  an  empirical  HR  Wiener 
filter  given  by 


kj- 


£<e<2V 

k~\ _ 

m 

Y.yf* 

k= 1 


m  = 


N_ 

V 


j  1,2, 


(7.16) 


where  N  is  the  signal  length,  j  is  the  level  of  the  DWT  decomposition,  and  J  represents  the 
maximum  decomposition  level.  This  empirical  DR  Wiener  filter  is  the  ratio  of  the  energy 
of  the  estimated  noise-free  signal  wavelet  coefficients  for  a  particular  scale  divided  by  the 
energy  of  the  noisy  signal  wavelet  coefficients  for  the  same  scale.  The  DR  Wiener  filter  is 
applied  to  a  DWT  of  the  original  noisy  signal  y<2>i  k,  formed  by  W3,  to  produce  the  wavelet 
coefficient  signal  estimate  .  An  inverse  DWT,  represented  by  W{1  is  applied  to  produce 
the  final  signal  estimate  s .  A  Daubechies-8  wavelet  is  used  for  the  DWT  in  block  Wj  while 
a  Daubechies-16  wavelet  is  used  for  the  DWT  in  blocks  W2  and  W3. 

Ghael  et  al  provide  the  following  explanation  for  the  improved  performance  of  their 
product  compared  to  traditional  thresholding  methods  [36].  The  goal  of  any  wavelet  based 
denoising  technique  is  to  determine  the  noise-free  signal  wavelet  coefficients  from  the  noisy 
signal  wavelet  coefficients.  The  noise-free  wavelet  coefficients  consist  of  a  number  of 
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trustworthy  coefficients  which  are  most  assuredly  due  to  the  signal  and  then  a  number  of  less 
trustworthy  coefficients  which  are  more  difficult  to  determine  due  to  the  presence  of  noise 
coefficients.  WienerShrink  uses  the  initial  thresholding  method  to  determine  the  trustworthy 
coefficients  given  by  The  untrustworthy  or  dubious  coefficients  are  then  predicted  by 
the  HR  Wiener  filter  formed  from  the  trustworthy  coefficients.  Ghael  et  al  indicate  that 
classical  wavelet  thresholding  methods  are  overly  conservative  in  thresholding  the 
coefficients  resulting  in  the  retention  of  only  the  trustworthy  coefficients. 

WienerShrink  is  applied  to  the  same  noisy  sinusoidal  signal  used  previously  with  the 
MSE  results  displayed  in  Figure  7.7  for  a  SNR  of  0  dB  and  Figure  7.8  for  a  SNR  of  10  dB. 
The  Heuristic  SURE  and  SURE  methods  of  thresholding  were  chosen  due  to  their  superior 
performance  over  minimax  and  universal  thresholding  methods  at  lower  SNR  levels.  The 
hard  threshold  option  outperforms  the  soft  threshold  option  in  some  areas,  especially  at  the 
lower  SNR  and  the  lowest  scale.  The  difference  is  not  substantial  allowing  either 
thresholding  option  to  be  utilized.  A  comparison  of  Figures  7.4  and  7.7  reveals  that  the 
minimum  MSE  obtained  by  both  the  DR  Wiener  wavelet  and  WienerShrink  at  a  particular 
scale  are  nearly  identical.  The  IIR  Wiener  wavelet  exhibits  superior  MSE  performance  in 
the  finer  scales  for  a  particular  decomposition  level.  For  example,  the  single-level  or  scale 
decomposition  for  HR  Wiener  wavelet  produces  an  average  MSE  value  of  0.37  in  the  lowest 
scale  (normalized  frequency  of  0.25  to  0.5)  while  the  best  WienerShrink  method  for  a  single- 
level  or  scale  decomposition  with  Heuristic  SURE  hard  thresholding  achieves  an  average 
MSE  value  of  only  0.6  at  the  same  scale. 
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If  the  filter  bank  can  be  adapted  by  changing  the  level  of  decomposition  (alternatively,  the 
number  of  scales)  to  confine  the  signal  of  interest  to  the  highest  scale  then  this  difference 
would  have  little  effect.  Otherwise,  the  HR  Wiener  wavelet  filter  may  be  preferable  to 
WienerShrink.  Thus  the  choice  becomes  a  signal  dependent  decision.  A  significant 
improvement  in  MSE  performance  for  Wienershrink  is  achieved  when  the  SNR  is  increased 
to  10  dB  (see  Figure  7.8).  This  difference  is  particularly  noticeable  for  denoising  methods 
which  employ  wavelet  thresholding. 

Comparing  wavelet  thresholding  methods  (Figure  5.9)  to  WienerShrink  (Figure  7.7) 
reveals  minimal  difference  between  the  two  when  applied  to  a  noisy  sinusoidal  signal. 
WienerShrink  is  slightly  better  in  terms  of  MSE  performance.  In  Chapter  VIII  signals  with 
more  interesting  characteristics  will  be  denoised  using  the  various  methods  and  the  results 
compared. 
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Heuristic  SURE  Soft 


Heuristic  SURE  Hard 


SURE  Soft 


SURE  Hard 


Figure  7.7:  WienerShrink  applied  to  noisy  sinusoidal  signal:  SNR  =  0  dB.  Four  level 
DWT  decomposition  using  Daubechies-8  to  form  pilot  signal  and  Daubechies-16  for  HR 
filtering.  Heuristic  SURE  and  SURE  thresholds  used  with  soft  and  hard  options. 
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Figure  7.8:  WienerShrink  applied  to  noisy  sinusoidal  signal:  SNR  =  10  dB.  Four  level 

DWT  decomposition  using  Daubechies-8  to  form  pilot  signal  and  Daubechies-16  for  DR 

filtering.  Heuristic  SURE  and  SURE  thresholds  used  with  soft  and  hard  options. 


VIH.  COMPARISON  OF  DENOISING  METHODS 

Various  Wiener  and  wavelet  filtering  methods  have  been  presented.  Thus  far  the 
ability  to  denoise  a  noisy  sinusoidal  signal  has  been  the  basis  for  comparison  between  the 
denoising  methods.  In  this  chapter  the  denoising  methods  are  applied  to  four  synthetic 
signals  which  are  commonly  used  as  benchmarks  for  comparison  because  of  their 
characteristics.  These  signals  were  originally  chosen  by  the  producers  of  Wavelab  at  the 
Stanford  University  Statistics  Department  [37].  The  noise-free  signals  are  shown  below  in 
Figure  8. 1  and  are  referred  to  as  Doppler,  HeaviSine,  Bumps,  and  Blocks. 
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Figure  8.1:  Noise-free  test  signals:  Doppler,  HeaviSine,  Bumps,  and  Blocks. 
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All  of  the  denoising  methods  discussed  in  this  thesis  are  applied  to  the  four  test 
signals  and  the  average  normalized  MSE  values  are  determined  for  twenty  trials.  The  results 
are  tabulated  in  the  Appendix.  Graphs  depicting  the  MSE  performance  of  the  various 
methods  are  displayed  as  Figures  8.2  to  8.8.  All  methods,  except  for  the  Wiener  filter,  are 
compared  in  Figures  8.9  and  8.10.  All  Figures  depicting  the  performance  of  Wiener  wavelet 
denoising  methods  show  two  different  assumptions  concerning  the  noise  estimation.  Figures 
8.4, 8.6,  and  8.9  assume  an  ideal  situation  in  which  the  exact  noise  sample  that  produced  the 
noisy  signal  is  supplied  to  the  Wiener  filter.  In  Figures  8.5,  8.7,  and  8.10,  an  independent 
white  Gaussian  noise  source  with  standard  deviation  equal  to  one  and  a  zero  mean  is 
provided  to  the  Wiener  filter.  The  second  option  is  realistic  and  better  approximates  the  true 
performance  of  these  methods. 

A  visual  comparison  of  four  denoising  methods  is  depicted  in  Figures  8.1 1  and  8. 12. 
Noise-free  and  noisy  (SNR=0  dB)  doppler  signals  of  length  1024  are  shown  in  Figure  8.11. 
The  noisy  doppler  signal  has  a  MSE  of  0.5774.  Denoised  doppler  results  are  shown  for 
wavelet  thresholding,  HR  Wiener  wavelet,  WaveShrink,  and  FIR  Wiener  wavelet  noise 
removal  techniques.  These  plots  provide  a  visual  interpretation  of  the  MSE  results.  Specific 
observations  regarding  the  denoising  methods  follow. 

A.  WIENER  FILTER 

The  Wiener  filter  did  not  denoise  these  signals  as  well  as  the  other  methods  except 
in  the  case  of  the  HeaviSine  signal.  Doppler,  Bumps,  and  Blocks  are  nonstationary  and 
therefore  are  not  effectively  filtered  by  a  standard  Wiener  filter.  HeaviSine  exhibits  some 
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degree  of  stationarity,  resulting  in  better  denoising  by  the  Wiener  filter.  A  short-time  Wiener 
filter  could  be  designed  with  an  optimal  window  size  and  filter  length  which  would  improve 
the  Wiener  filter  MSE  performance  significantly.  This  type  of  filter  is  utilized  in  denoising 
nonstationary  signals  in  the  next  chapter. 

B.  WAVELET  THRESHOLDING 

Wavelet  thresholding  or  “Wave  Shrink”  performs  as  well  as  or  better  than  all  the 
other  methods  at  SNR  levels  from  -5  to  20  dB.  In  general  the  SURE  and  Hybrid  thresholding 
methods  are  the  best  of  the  four  thresholding  methods.  The  appendix  tables  show  results  for 
both  soft  and  hard  threshold  options,  although  only  the  soft  option  was  used  in  the  plots.  In 
some  cases  better  MSE  results  are  obtained  using  the  hard  threshold,  however,  the  soft 
threshold  provided  better  overall  results.  The  signals  in  this  study  have  all  been  scaled  in 
amplitude  such  that  the  desired  SNR  is  achieved  while  maintaining  the  noise  level  at  a 
standard  deviation  of  one.  Recall  the  method  for  determining  the  noise  standard  deviation 
in  actual  signals  required  determining  the  mean  absolute  deviation  from  the  wavelet 
coefficients  at  the  lowest  or  finest  scale.  In  cases  where  some  of  the  signal  energy  resides 
in  this  finest  scale,  the  performance  of  the  wavelet  threshold  denoising  method  will  be 
degraded. 

C.  FIR  WIENER  WAVELET  FILTER 

The  FIR  Wiener  wavelet  filter  is  the  best  denoising  method  in  terms  of  MSE 
performance  in  the  ideal  noise  conditions  described  above  and  shown  in  Figures  8.4  and  8.9. 
Figures  8.5  and  8.10  show  that  this  method  is  affected  the  most  when  provided  with  an 
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independent  noise  source.  The  MSE  more  than  doubles  for  some  of  the  test  signals.  An 
accurate  noise  estimate  is  critical  to  the  performance  of  the  FIR  Wiener  wavelet  denoising 
method.  Transient  affects  also  must  be  considered  when  using  a  FIR  filter.  These 
shortcomings  make  this  method  the  least  preferred. 

D.  HR  WIENER  WAVELET  FILTER 

The  DR  Wiener  wavelet  filter  exhibits  a  MSE  performance  that  is  comparable  to 
wavelet  thresholding  and  WienerShrink.  Although  the  MSE  error  performance  is  degraded 
when  assuming  an  independent  noise  estimate  as  seen  by  comparing  Figures  8.6  and  8.7,  it 
still  provides  substantial  denoising  capability.  Furthermore,  it  is  the  easiest  to  implement, 
requires  the  least  computational  effort,  and  has  no  transient  effects. 

E.  WIENERSHRINK 

WienerShrink  outperformed  all  other  methods  but  not  by  a  substantial  amount.  It  is 
likely  that  Ghael  et  al,  the  developers  of  this  technique,  may  have  applied  an  unpublished 
statistical  thresholding  method  in  their  original  work  which  was  not  replicated  in  this  thesis. 
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Figure  8.2:  MSE  vs.  SNR  performance  for  Wiener  Filter  (Filter  Order  8,12,16,20)  on  four 
test  signals.  Twenty  trial  average. 


*  SURE 
+ 


HeaviSine 


co  0.15 


o  Mini 


co  0.15 


5  10  15 

0 

-5  0 

5  10  1 

5 

SNR  (dB) 

SNR  (dB) 

Bumps 

Blocks 

co  0.15 


UJ  0.08 


.0 


SNR  (dB)  SNR  (dB) 

Figure  8.3:  MSE  vs.  SNR  performance  for  wavelet  threshold  denoising  of  four  test 
signals.  Threshold  types  shown  include  SURE,  Heuristic  SURE,  Minimax,  and 
Universal.  Five  level  or  scale  decomposition  using  Daubechies~20  wavelet.  Twenty  trial 
average. 
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Figure  8.4:  MSE  vs.  SNR  performance  for  FIR  Wiener  wavelet  denoising  of  four  test 
signals.  Wiener  filter  orders  of  8  and  16  shown.  Five  level  or  scale  decomposition  using 
Daubechies-20  wavelet.  Twenty  trial  average. 
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Figure  8.5:  MSE  vs.  SNR  performance  for  FIR  Wiener  wavelet  denoising  of  four  test 
signals.  Wiener  filter  provided  independent  noise  source.  Filter  orders  of  8  and  16 
shown.  Five  level  or  scale  decomposition  using  Daubechies-20  wavelet.  Twenty  trial 
average. 
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Figure  8.6:  MSE  vs.  SNR  performance  for  OR  Wiener  wavelet  denoising  of  four  test 
signals.  Three,  four,  and  five  level  or  scale  decomposition  shown  using  Daubechies-20 
wavelet.  Twenty  trial  average. 
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Figure  8.7:  MSE  vs.  SNR  performance  for  HR  Wiener  wavelet  denoising  of  four  test 
signals.  Wiener  filter  provided  with  independent  noise  source.  Three,  four,  and  five  level 
or  scale  decomposition  shown  using  Daubechies-20  wavelet.  Twenty  trial  average. 
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Figure  8.8:  MSE  vs.  SNR  performance  for  WienerShrink  denoising  of  four  test  signals 
SURE,  Heuristic  SURE,  Minimax,  and  Universal  thresholding  used  with  soft  option. 
Five  level  or  scale  decomposition  with  Daubechies-20  wavelet.  Twenty  trial  average. 
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Figure  8.10:  MSE  vs.  SNR  performance  comparison  of  wavelet  thresholding, 
WienerShrink,  FIR  Wiener  wavelet,  and  IIR  Wiener  wavelet  denoising  of  four  test 
signals.  Wiener  filters  provided  with  independent  noise  source. 


Noise-Free  Doppler 


Figure  8.11:  Noise-free  and  noisy  doppler  signal  of  length  1024  data  points.  SNR=0  dB. 
MSE  of  noisy  doppler  is  0.5774. 
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WaveShrink:  MSE=0.1021 


FIR  Wien/Wave:  MSE=0.1449 


Figure  8.12:  Comparison  of  four  denoising  techniques  applied  to  noisy  doppler  with 
SNR=0  dB.  MSE  values  are  shown  for  each  method.  Independent  noise  sources 
provided  to  Wiener  filters.  Daubechies-20  wavelet  used  in  five-level  or  scale 
decompositions.  SURE  (soft)  threshold  used  for  wavelet  thresholding. 
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IX.  WAVELET  PACKETS 


Wavelet  analysis  has  been  shown  to  be  particularly  effective  at  removing  noise  from 
low  frequency  signals  such  as  narrowband  tonals  which  might  occur  in  the  ocean 
environment.  The  constant-Q  filters  provide  a  logarithmic  frequency  resolution  with  narrow 
bandwidths  at  low  frequencies  and  wide  bandwidths  at  high  frequencies.  This  high 
frequency  resolution  at  low  frequencies  is  produced  by  iterating  the  lowpass  scaling  branch 
of  the  Mallat  algorithm  tree.  In  addition  to  low  frequency  tonals,  ocean  acoustic  signals  often 
include  high  frequency  signals  which  may  be  of  interest  as  well.  These  signals  are  usually 
short  duration  transients  or  acoustic  energy  pulses.  Wavelet  based  denoising  techniques  are 
not  optimal  for  removing  noise  from  this  class  of  signals. 

A  richer  signal  analysis  tool,  termed  wavelet  packets,  was  introduced  by  Ronald 
Coifman  [38]  to  provide  high  resolution  decomposition  at  high  frequencies.  This  is 
accomplished  by  iterating  the  highpass  wavelet  branch  of  the  Mallat  algorithm  tree.  The 
detail  coefficient  vector  is  thus  decomposed  into  two  parts  by  splitting,  filtering  and 
decimating  in  the  same  manner  as  the  approximation  coefficient  vector.  The  full  binary  tree 
is  pictured  in  Figure  9.1  where  the  H  and  L  blocks  are  the  high  and  low-pass  filters 
determined  by  the  wavelet  and  scaling  functions.  The  numerous  signal  expansions  that  are 
possible  with  wavelet  packets  come  at  a  cost  in  computational  complexity  of  0(N  log l0(N)) 
compared  to  0(N)  for  the  wavelet  transform. 
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Figure  9.1:  Wavelet  packet  tree  (Used  with  permission  from  Joshua  Altmann,  [14]). 


A.  WAVELET  PACKET  THEORY 

Wavelet  packets  are  formed  by  modifying  (4.16 )  and  (4.17 ),  which  determine  the 
DWT  filter  coefficients  in  terms  of  the  scaling  and  wavelet  functions.  A  third  index  is 
required  to  describe  the  wavelet  packet  in  terms  of  its  position  within  the  wavelet  packet  tree. 
This  index  describes  the  wavelet  packet  in  terms  of  its  frequency  resolution.  The  wavelet 
packet  equations  are  given  by 


=  n  e  Z 

n 

typ(0  -  ]Tg(n)\/2<t>p(2t-«),  n  6  Z. 


(9.1) 


The  wavelet  packets  for  a  three-level  wavelet  packet  decomposition  using  the  familiar 
Daubechies-2  wavelet  are  shown  in  Figure  9.2.  These  wavelet  packets  correspond  to  the 
Figure  9.1  tree  with  the  notation  WJ  p  with  j  =  3  denoting  the  scale  parameter  and  p  =  0  to 
7  the  frequency  parameter.  Notice  the  first  two  wavelet  packets  correspond  to  the 
Daubechies-2  filters  (Figure  4.1 1)  developed  from  the  scaling  and  wavelet  functions  for  the 
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Figure  9.2:  Daubechies-2  wavelet  packets  for  a  three-level  or  scale  decomposition.  WhV 
for  j  =  3  and  p  =  0  to  7.  Supported  on  the  interval  [0, 31- 

wavelet  transform.  Thus  the  wavelet  transform  is  a  subset  of  wavelet  packets. 

The  wavelet  packet  decomposition  of  a  signal  results  in  a  binary  tree  composed  of 

waveforms  known  as  “atoms.”  This  collection  of  atoms  forms  an  overcomplete 

“dictionary”  or  library  from  which  a  signal  of  length  N  can  be  decomposed  and  reconstructed 

in  at  most  2N  different  ways.  The  Best  Basis  algorithm,  developed  by  R.  Coifman  and  V. 

Wickerhauser  [38],  and  used  in  MATLAB  [22],  seeks  to  minimize  an  additive  cost  function 

in  order  to  determine  the  best  basis  for  representing  the  signal.  The  cost  function  chosen  is 
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entropy  given  by 


E(s)  =  sf  log (sf)  (9.3) 

i 

where  s  is  the  signal  and  st  the  coefficients  of  s  in  an  orthonormal  basis.  Best  Basis  is  an 
efficient  algorithm  which  computes  the  entropy  at  each  node  in  the  binary  tree  and  then  finds 
the  best  wavelet  packet  representation  of  the  signal  based  on  minimizing  this  quantity. 
Reference  [34]  provides  additional  information. 

B.  WAVELET  PACKETS  APPLIED  TO  TEST  SIGNAL 

The  test  signal  chosen  for  analysis  is  a  series  of  three  high  frequency  transient  pulses 
of  decreasing  amplitude.  This  test  signal,  which  will  be  referred  to  as  transients,  was  used 
by  Barsanti  in  reference  [8].  The  signal  is  produced  from  the  following  truncated 
exponentially  decaying  sinusoids: 

51  =  sin(2n:£7/4)  *exp(  -£7/200)  kl  =  1,2,..., 256 

52  =  1/2  sin(27t£2/4)*exp( -£2/200)  £2  =  1,2,...,200  (9.2) 

53  =  1/3  sin(27t£3/4)  *exp( -£3/200)  £3  =  1,2,..., 128  . 

This  is  a  realistic  acoustic  signal  which  is  artificially  produced  to  allow  MSE  comparisons 

at  differing  SNR  levels.  The  noisy  version  of  transients  is  produced  by  adding  white 
Gaussian  noise  with  a  zero  mean.  The  standard  deviation  of  the  noise  remains  constant  at 
a  value  of  one  while  the  noise-free  signal  is  scaled  to  obtain  the  desired  SNR  level.  Noise- 
free  transients  and  noisy  transients  with  an  SNR  of  0  dB  are  shown  in  Figure  9.3.  The  first 
transient  is  also  pictured  expanded  in  the  noise-free  and  noisy  form  so  that  the  details  may 
be  better  observed. 
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Noise-free  Transients 


Single  Noise-free  Transient 


Noisy  Transients:  SNR  =  0  dB  Single  Noisy  Transient 
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Figure  9.3:  Noise-free  and  noisy  transients  (SNR  =  0  dB).  Expanded  view  of  first 
transient  without  and  with  noise  (SNR  =  0  dB). 


Seven  different  methods  of  denoising  are  applied  to  the  transients.  The  results  for 
the  short-time  Wiener  filter  is  shown  in  Figure  9.4.  The  other  six  methods  are  compared  in 
Figures  9.5  through  9.10  for  SNR  levels  of  5, 0,  and  -5  dB.  Specific  comments  concerning 
each  method  follow. 
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1. 


Short-time  Wiener  Filter 


The  short-time  Wiener  filter  of  length  twenty  is  applied  to  the  16384  point  transient 
signal  in  segments  of  256  points.  The  noise  source  provided  to  the  Wiener  filter  is  an 
independent  noise  source  from  that  used  to  produce  the  noisy  transients  signal.  It  is  white 
Gaussian  noise  with  standard  deviation  equal  to  one  and  a  zero  mean.  Windowing  the  data 
allows  the  Wiener  filter  to  completely  remove  the  noise  in  noise-only  segments.  This  does 
not  occur  when  the  Wiener  filter  is  applied  to  the  entire  signal  in  one  block  due  to  the  non¬ 
stationary  nature  of  the  data. 

2.  Wavelet  Thresholding 

This  signal  does  not  meet  the  low  frequency  assumption  of  wavelet  denoising  and 
therefore  the  MSE  is  greater  than  that  achieved  by  wavelet  packet  thresholding.  When  using 
wavelet  thresholding  for  denoising,  the  approximation  coefficients  are  not  usually 
thresholded  since  the  assumption  is  that  these  coefficients  contain  the  signal  for  the  most 
part.  The  algorithm  was  modified  to  support  approximation  coefficient  thresholding  which 
removed  additional  noise  and  reduced  the  MSE. 

3.  FIR  Wiener  Wavelet  Filter 

The  FIR  Wiener  wavelet  filter  does  not  perform  well  for  this  signal,  again  due  to  the 
high  frequency  characteristics  of  the  signal.  The  MSE  is  higher  than  what  could  be  optimally 
achieved  due  to  the  residual  noise  present  in  the  noise-only  segments.  This  is  true  for  all  of 
the  Wiener/Wavelet  methods.  This  is  partially  remedied  by  implementing  a  windowed  or 
segmented  filtering  algorithm  that  uses  a  triangular  window  with  fifty-  percent  overlap.  A 
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window  length  of  4096  points  with  a  Wiener  filter  length  of  four  and  decomposition  level 
of  two  achieved  the  best  results  which  are  listed  in  the  table  at  the  end  of  this  section.  The 
short-time  FIR  Wiener  wavelet  filter  approaches  the  MSE  performance  of  the  FIR  Wiener 
Wavelet  packet  filter. 

4.  HR  Wiener  Wavelet  Filter 

The  observations  made  for  the  FIR  Wiener  wavelet  filter  apply  to  this  method  as  well. 
The  DR  Wiener/Wavelet  filter  experiences  the  most  severe  degradation  when  the  signal  to 
be  denoised  doesn’t  meet  the  low  frequency  assumption.  Significant  improvement  is 
achieved  by  implementing  the  same  triangular  window  short-time  filtering  algorithm 
described  above  for  the  short-time  FIR  Wiener  wavelet  filter.  A  window  length  of  256 
points  was  used  to  produce  the  results  shown  in  Figure  9.11.  Additional  improvement  in 
terms  of  MSE  is  realized  by  applying  a  threshold  to  the  HR  Wiener  wavelet  filter 
coefficients.  If  the  filter  coefficient  magnitude  is  below  some  threshold  value,  in  this  case 
0.4  was  chosen  by  trial  and  error,  then  the  coefficient  is  set  to  zero.  The  resulting 
improvement  is  seen  in  column  four  of  Figure  9. 1 1. 

5.  Wavelet  Packet  Thresholding 

Wavelet  Packet  thresholding  exhibits  significant  MSE  improvement  over  wavelet 
thresholding.  Greater  frequency  resolution  at  high  frequencies  allows  this  method  to  produce 
transient  pulses  with  better  fidelity.  Wavelet  packet  thresholding  uses  the  SURE  threshold 
described  previously  in  Chapter  IV.  The  approximation  coefficients,  represented  by  the 
“atom”  in  the  lower  left  comer  of  the  binary  tree,  are  thresholded  as  well,  leading  to  a  noise- 
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free  result  in  the  segments  between  the  pulses  which  contained  only  noise. 

6.  FIR  Wiener  Wavelet  Packet  Filter 

In  this  method  a  separate  FIR  Wiener  filter  is  applied  to  each  atom  in  the  terminal  or 
lowest  nodes  of  the  binary  tree.  No  Best  Basis  optimization  is  performed.  Some  MSE 
improvement  is  realized  when  comparing  this  method  to  the  FIR  Wiener  wavelet  method  but, 
the  short-time  Wiener  filter  and  wavelet  packet  thresholding  still  achieve  better  results. 

7.  IIR  Wiener  Wavelet  Packet  Filter 

In  this  method  an  DR  Wiener  Wavelet  filter  is  applied  to  each  atom  at  the  terminal 
nodes.  The  IIR  Wiener  Wavelet  Packet  filter  displays  some  improvement  over  the  DR 
Wiener  wavelet  filter  but  it  is  the  poorest  performer  of  all  the  wavelet  packet  methods.  It 
would  likely  achieve  better  results  in  some  type  of  windowed  or  segmented  implementation. 

8.  Summary 

The  table  following  shows  a  comparison  of  all  the  methods  in  increasing  MSE  order. 
An  optimized  ST  Wiener  filter  still  outperforms  all  other  methods  with  the  wavelet  packet 
threshold  method  and  the  ST  DR  Wiener/wavelet  method  following  close  behind.  Methods 
which  do  not  utilize  wavelet  packets  or  a  short-time  filter  are  unable  to  effectively  denoise 
this  high  frequency  transient  signal.  For  this  particular  signal,  a  short-time  implementation 
of  wavelet  packet  based  denoising  methods  produced  no  further  improvement  and  was 
computationally  cumbersome. 
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|  Method 

-5  dB 

0  dB 

5  dB 

ST  Wiener 

0.1468 

0.0684 

0.0363 

Wavelet  Packet  Threshold 

0.2888 

0.1231 

0.0439 

ST  HR  Wiener/wavelet 

0.3751 

0.1562 

0.0618 

Wavelet  Packet  FIR  Wiener 

0.4102 

0.1884 

0.0969 

ST  FIR  Wiener/wavelet 

0.5485 

0.2490 

0.1115 

FIR  Wiener/wavelet 

0.6944 

0.3558 

0.1768 

Wavelet  Packet  HR  Wiener 

0.7055 

0.3029 

0.1283 

Wavelet  Threshold 

0.8283 

0.3832 

0.1399 

DR  Wiener/wavelet 

1.264 

0.8555 

0.4788 
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Short-time  Wiener  (SNR=-5dB):  MSE=.1468 
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Short-time  Wiener  (SNR=0dB):  MSE=.0684 


Short-time  Wiener  (SNR=5dB):  MSE=.0363 


Figure  9.4:  Denoising  transients  using  short-time  Wiener  filter  of  length  twenty. 
Window  size:  256.  SNR  levels  of  -5, 0,  5  dB. 
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Figure  9.5:  Denoising  transients  using  wavelet  thresholding,  FIR  Wiener/Wavelet,  and 
HR  Wiener/Wavelet  methods.  SNR  =  5  dB. 
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Figure  9.6:  Denoising  transients  using  wavelet  packet  thresholding,  FIR 
Wiener/Wavelet  packet,  and  HR  Wiener/Wavelet  packet  methods.  SNR  =  5  dB. 
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Figure  9.7:  Denoising  transients  using  wavelet  thresholding,  FIR  Wiener/Wavelet, 
and  DR  Wiener/Wavelet  methods.  SNR  =  0  db. 
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Noisy  Transients:  SNR=OdB 
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WP  Thresh:  MSE=.1231 
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WP  HR  Wien/Wave:  MSE=.3029 


Figure  9.8:  Denoising  transients  using  wavelet  packet  thresholding,  FIR 
Wiener/Wavelet  packet,  and  HR  Wiener/Wavelet  packet  methods.  SNR  =  0  dB. 
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FIR  Wien/Wave:  MSE=0.6944 


HR  Wien/Wave:  MSE=1.2640 


Figure  9.9:  Denoising  transients  using  wavelet  thresholding,  FIR  Wiener/Wavelet, 
and  HR  Wiener/Wavelet  methods.  SNR  =-  5  dB. 
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Figure  9.10:  Denoising  transients  using  wavelet  packet  thresholding,  FIR 
Wiener/Wavelet  packet,  and  HR  Wiener/Wavelet  packet  methods.  SNR  =  -5  dB. 
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Figure  9.1 1:  Short-time  HR  Wiener  wavelet  filter  applied  to  transients  with  SNR  levels  of 
5,  0,  and  -5  dB.  1st  column:  noisy  transients.  2nd  column:  denoised  transients.  3rd 
column:  denoised  transients  using  ST  HR  Wiener  wavelet  filter  with  hard  threshold. 
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X.  CONCLUSIONS 


This  study  was  undertaken  to  determine  the  feasibility  of  combining  Wiener  filter  and 
wavelet  based  denoising  techniques.  The  two  methods  have  been  used  independently  for 
years  to  denoise  many  types  of  signals.  Their  development  and  performance  were 
demonstrated  in  Chapters  HI  through  V.  Before  any  attempt  was  made  to  combine  the 
methods,  legitimate  issues  arose  concerning  the  aliasing  effects  produced  by  inserting  a 
Wiener  filter  into  the  DWT  filterbank  and  how  this  might  affect  the  ability  to  perfectly 
reconstruct  the  denoised  signal.  This  concern  led  to  the  development  of  an  alias  cancellation 
filter  in  Chapter  VI.  It  was  hoped  that  this  would  cancel  the  deleterious  affects  produced  by 
aliasing.  Unfortunately,  the  alias  cancellation  filter  only  produced  greater  distortion  and  was 
not  found  to  be  of  any  benefit. 

Though  the  original  intent  was  to  combine  a  FIR  Wiener  filter  with  wavelet  analysis, 
the  possibility  of  an  HR  Wiener  wavelet  filter  was  motivated  by  the  WienerShrink  method 
of  reference  [36].  The  simplicity  of  this  method  combined  with  its  advertised  performance 
led  to  the  DR  Wiener  wavelet  filter  methods  developed  in  Chapter  VII. 

In  the  earlier  chapters  the  Wiener  filter,  wavelet  thresholding,  the  FIR  Wiener 
wavelt  filter  and  DR  Wiener  wavelet  filters  were  developed  and  applied  to  a  stationary  noisy 
sinusoidal  signal.  Initial  results  indicated  that  the  HR  and  FIR  Wiener  wavelet  filters  might 
actually  outperform,  in  a  MSE  sense,  the  more  traditional  Wiener  filter  and  wavelet 
thresholding  methods.  These  results  proved  that  the  Wiener  wavelet  combination  could 
indeed  denoise  a  signal.  However,  removing  the  noise  from  a  single  sinusoidal  waveform 
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was  not  deemed  to  be  a  comprehensive  test  of  the  denoising  capabilities  of  these  techniques. 
Therefore,  the  robustness  of  the  various  methods  was  measured  by  applying  them  to  four  test 
signals.  In  these  results,  the  DR  Wiener  wavelet  methods  compared  favorable  with  wavelet 
thresholding.  The  HR  Wiener  wavelet  filter  results  which  seemed  promising  in  Chapter  VI 
were  rather  disappointing  compared  to  the  other  methods  with  MSE  values  as  much  as  twice 
those  of  the  other  methods. 

Since  many  ocean  acoustic  signals  of  interest  are  high  frequency  transients,  wavelet 
packet  methods  were  developed  in  Chapter  IX  to  provide  greater  frequency  resolution  at 
high  frequencies.  The  HR  and  HR  Wiener  filtering  techniques  were  applied  to  the  terminal 
nodes  of  the  wavelet  packet  decomposition  tree  providing  an  enhanced  capability  in 
denoising  high  frequency  signals  compared  to  the  previously  developed  Wiener  wavelet 
methods.  The  various  denoising  methods  were  applied  to  synthetically  generated  noisy 
transients.  The  results  were  similar  in  that  the  traditional  methods  of  short-time  Wiener 
filtering  and  wavelet  packet  thresholding  outperformed  all  other  methods.  However,  an  HR 
Wiener  wavelet  filter  used  in  a  short-time  implementation  performed  rather  well  with  MSE 
values  somewhat  greater  than  the  traditional  methods. 

This  study  has  proven  the  feasibility  of  combining  Wiener  and  wavelet  based 
techniques  into  a  single  filter.  In  general,  it  does  not  outperform  either  of  the  methods  from 
which  it  is  derived.  The  basis  for  comparison  between  the  various  methods  is  not  by  any 
means  clear  cut.  Each  method  has  numerous  parameters  which  can  be  adjusted  to  optimize 
its  denoising  performance.  To  complicate  matters  further,  many  of  these  parameters  are 
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signal  dependent.  Although  the  conclusions  drawn  are  well  supported,  specific  results  in 
denoising  a  particular  signal  may  vary  depending  on  a  variety  of  parameters  such  as  Wiener 
filter  length,  wavelet  filter  length,  window  length,  wavelet  type,  threshold  type,  hard  or  soft 
threshold  option,  thresholding  of  approximation  coefficients,  and  the  list  continues.  The 
variability  makes  this  a  fascinating  topic  which  often  leads  to  more  questions  than  answers. 
No  doubt  this  subject  will  stimulate  interest  for  years  to  come.  Future  studies  might  involve 
the  application  of  a  median  filter  to  the  wavelet  coefficients  or  some  type  of  hybrid  Wiener 
wavelet  threshold  method  in  which  a  minimum  filter  coefficient  value  is  set  based  on 
statistical  analysis. 
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APPENDIX 


This  appendix  lists  modified  MSE  results  as  defined  by  (3.1 1).  The  MSE  results  are 
tabulated  for  the  four  test  signals  Doppler,  HeaviSine,  Bumps,  and  Blocks  which  are  pictured 
in  Figure  8.1.  Variable  parameters  include:  SNR,  Wiener  filter  order,  level  of 
decomposition,  wavelet  threshold  type,  and  wavelet  threshold  option  (s  for  soft,  h  for  hard). 
Actual  and  independent  noise  sample  refers  to  the  noise  provided  to  the  Wiener  filter.  In  the 
actual  case  the  noise  used  to  produce  the  noisy  signal  provided,  while  the  independent  case 
refers  to  an  independent  sample  of  white  Gaussian  noise  with  a  zero  mean  and  variance  as 
required  to  meet  the  SNR  level  desired. 
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